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 "InterpolationCU.txx"
39 #include "InterpolationCC.txx"
41 using namespace MEDCoupling;
43 MEDCouplingRemapper::MEDCouplingRemapper():_src_ft(0),_target_ft(0),_interp_matrix_pol(IK_ONLY_PREFERED),_nature_of_deno(NoNature),_time_deno_update(0)
47 MEDCouplingRemapper::~MEDCouplingRemapper()
53 * This method is the second step of the remapping process. The remapping process works in three phases :
55 * - Set remapping options appropriately
56 * - The computation of remapping matrix
57 * - Apply the matrix vector multiply to obtain the result of the remapping
59 * 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.
60 * 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).
61 * So this method is less precise but a more user friendly way to compute a remapping matrix.
63 * \param [in] srcMesh the source mesh
64 * \param [in] targetMesh the target mesh
65 * \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.
66 * 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).
68 * \sa MEDCouplingRemapper::prepareEx
70 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
72 if(!srcMesh || !targetMesh)
73 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepare : presence of NULL input pointer !");
74 std::string srcMethod,targetMethod;
75 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
76 MCAuto<MEDCouplingFieldTemplate> src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
77 src->setMesh(srcMesh);
78 MCAuto<MEDCouplingFieldTemplate> target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
79 target->setMesh(targetMesh);
80 return prepareEx(src,target);
84 * This method is the generalization of MEDCouplingRemapper::prepare. Indeed, MEDCouplingFieldTemplate instances gives all required information to compute the remapping matrix.
85 * This method must be used instead of MEDCouplingRemapper::prepare if Gauss point to Gauss point must be applied.
87 * \param [in] src is the field template source side.
88 * \param [in] target is the field template target side.
90 * \sa MEDCouplingRemapper::prepare
92 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
95 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL input pointer !");
96 if(!src->getMesh() || !target->getMesh())
97 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
99 _src_ft=const_cast<MEDCouplingFieldTemplate *>(src); _src_ft->incrRef();
100 _target_ft=const_cast<MEDCouplingFieldTemplate *>(target); _target_ft->incrRef();
101 if(isInterpKernelOnlyOrNotOnly())
102 return prepareInterpKernelOnly();
104 return prepareNotInterpKernelOnly();
107 int MEDCouplingRemapper::prepareInterpKernelOnly()
109 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
117 // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED = 10,
118 // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED = 11,
120 // } MEDCouplingMeshType;
122 switch(meshInterpType)
124 case 90: // UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
125 case 91: // UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
126 case 165: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
127 case 181: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - UNSTRUCTURED
128 case 170: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
129 case 171: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
130 case 186: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
131 case 187: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
132 case 85: // UNSTRUCTURED - UNSTRUCTURED
133 return prepareInterpKernelOnlyUU();
134 case 167: // SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
135 case 183: // SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED - CARTESIAN
136 case 87: // UNSTRUCTURED - CARTESIAN
137 return prepareInterpKernelOnlyUC();
138 case 122: // CARTESIAN - SINGLE_STATIC_GEO_TYPE_UNSTRUCTURED
139 case 123: // CARTESIAN - SINGLE_DYNAMIC_GEO_TYPE_UNSTRUCTURED
140 case 117: // CARTESIAN - UNSTRUCTURED
141 return prepareInterpKernelOnlyCU();
142 case 119: // CARTESIAN - CARTESIAN
143 return prepareInterpKernelOnlyCC();
144 case 136: // EXTRUDED - EXTRUDED
145 return prepareInterpKernelOnlyEE();
147 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
151 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
153 std::string srcm,trgm,method;
154 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
155 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
158 return prepareNotInterpKernelOnlyGaussGauss();
161 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
162 throw INTERP_KERNEL::Exception(oss.str().c_str());
168 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
169 * 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.
171 * \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.
172 * \param [in/out] targetField the destination field with the allocated array in which all tuples will be overwritten.
173 * \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.
174 * 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
175 * 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
176 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
180 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
182 if(!srcField || !targetField)
183 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transfer : input field must be both not NULL !");
184 transferUnderground(srcField,targetField,true,dftValue);
188 * This method is equivalent to MEDCoupling::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
189 * 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
191 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
193 * \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.
194 * \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.
196 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
198 if(!srcField || !targetField)
199 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : input field must be both not NULL !");
200 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
203 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
205 if(!srcField || !targetField)
206 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::reverseTransfer : input fields must be both not NULL !");
208 targetField->checkConsistencyLight();
209 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
210 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
211 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
212 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
213 if(srcField->getNature()!=targetField->getNature())
214 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
215 if(targetField->getNumberOfTuplesExpected()!=_target_ft->getNumberOfTuplesExpected())
217 std::ostringstream oss;
218 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();
219 oss << " ! It appears that the target support is not the same between the prepare and the transfer !";
220 throw INTERP_KERNEL::Exception(oss.str().c_str());
222 DataArrayDouble *array(srcField->getArray());
223 int trgNbOfCompo=targetField->getNumberOfComponents();
226 srcField->checkConsistencyLight();
227 if(trgNbOfCompo!=srcField->getNumberOfTuplesExpected())
228 throw INTERP_KERNEL::Exception("Number of components mismatch !");
232 MCAuto<DataArrayDouble > tmp(DataArrayDouble::New());
233 tmp->alloc(srcField->getNumberOfTuplesExpected(),trgNbOfCompo);
234 srcField->setArray(tmp);
236 computeDeno(srcField->getNature(),srcField,targetField);
237 double *resPointer(srcField->getArray()->getPointer());
238 const double *inputPointer=targetField->getArray()->getConstPointer();
239 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
243 * This method performs the operation source to target using matrix computed in MEDCoupling::MEDCouplingRemapper::prepare method.
244 * If mesh of \b srcField does not match exactly those given into \ref MEDCoupling::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
246 * \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.
247 * \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.
248 * 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
249 * 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
250 * be intercepted by a source cell on a different process. In this case 0. assigned to \a dftValue is more appropriate.
251 * \return destination field to be deallocated by the caller.
255 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
259 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input srcField is NULL !");
260 srcField->checkConsistencyLight();
261 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
262 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
263 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
264 ret->setNature(srcField->getNature());
265 transfer(srcField,ret,dftValue);
266 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
270 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
273 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferField : input targetField is NULL !");
274 targetField->checkConsistencyLight();
276 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
277 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
278 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
279 ret->setNature(targetField->getNature());
280 reverseTransfer(ret,targetField,dftValue);
281 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
286 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
287 * is here only for automatic CORBA generators.
289 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
291 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
295 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
296 * is here only for automatic CORBA generators.
298 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
300 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
304 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
305 * is here only for automatic CORBA generators.
307 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
309 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
313 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or prefered.
314 * If interpolation matrix policy is :
316 * - 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
317 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
318 * If not, the \b not only INTERP_KERNEL method will be attempt.
320 * - 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
321 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
322 * If not, the INTERP_KERNEL only method will be attempt.
324 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
326 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
328 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
330 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
332 return _interp_matrix_pol;
336 * 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
337 * 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
338 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
340 * If interpolation matrix policy is :
342 * - 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
343 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
344 * If not, the \b not only INTERP_KERNEL method will be attempt.
346 * - 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
347 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
348 * If not, the INTERP_KERNEL only method will be attempt.
350 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
352 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
354 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c MEDCoupling::InterpolationMatrixPolicy
355 * for automatic generation of CORBA component.
357 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
359 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
361 switch(newInterpMatPol)
364 _interp_matrix_pol=IK_ONLY_PREFERED;
367 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
370 _interp_matrix_pol=IK_ONLY_FORCED;
373 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
376 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 !");
380 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
382 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
383 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
384 std::string srcMeth,trgMeth;
385 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
386 const int srcMeshDim=src_mesh->getMeshDimension();
389 srcSpaceDim=src_mesh->getSpaceDimension();
390 const int trgMeshDim=target_mesh->getMeshDimension();
393 trgSpaceDim=target_mesh->getSpaceDimension();
394 if(trgSpaceDim!=srcSpaceDim)
395 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
396 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
398 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
400 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
401 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
402 INTERP_KERNEL::Interpolation1D interpolation(*this);
403 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
405 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
407 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
408 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
409 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
410 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
412 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
414 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
415 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
416 INTERP_KERNEL::Interpolation2D interpolation(*this);
417 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
419 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
421 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
422 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
423 INTERP_KERNEL::Interpolation3D interpolation(*this);
424 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
426 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
428 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
429 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
430 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
431 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
433 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
435 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
436 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
437 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
438 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
439 INTERP_KERNEL::Interpolation3D interpolation(*this);
440 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
442 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
444 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
445 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
446 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
447 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
448 INTERP_KERNEL::Interpolation3D interpolation(*this);
449 std::vector<std::map<int,double> > matrixTmp;
450 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
451 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
452 ReverseMatrix(matrixTmp,nbCols,_matrix);
453 nbCols=matrixTmp.size();
455 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
457 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
459 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
460 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
461 INTERP_KERNEL::Interpolation2D interpolation(*this);
462 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
466 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
467 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
468 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
469 std::vector<std::map<int,double> > matrixTmp;
470 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
471 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
472 ReverseMatrix(matrixTmp,nbCols,_matrix);
473 nbCols=matrixTmp.size();
474 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
475 if(!duplicateFaces.empty())
477 std::ostringstream oss; oss << "An unexpected situation happened ! For the following 1D Cells are part of edges shared by 2D cells :\n";
478 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
480 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
481 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
487 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
489 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
491 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
492 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
493 INTERP_KERNEL::Interpolation2D interpolation(*this);
494 std::vector<std::map<int,double> > matrixTmp;
495 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
496 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
497 ReverseMatrix(matrixTmp,nbCols,_matrix);
498 nbCols=matrixTmp.size();
502 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
503 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
504 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
505 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
506 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
507 if(!duplicateFaces.empty())
509 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
510 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
512 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
513 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
519 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
521 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
522 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
523 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
524 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
525 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
526 if(!duplicateFaces.empty())
528 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
529 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
531 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
532 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
537 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
539 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
540 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
541 INTERP_KERNEL::Interpolation2D3D interpolation(*this);
542 std::vector<std::map<int,double> > matrixTmp;
543 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
544 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
545 ReverseMatrix(matrixTmp,nbCols,_matrix);
546 nbCols=matrixTmp.size();
547 INTERP_KERNEL::Interpolation2D3D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
548 if(!duplicateFaces.empty())
550 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
551 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
553 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
554 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
559 else if(trgMeshDim==-1)
561 if(srcMeshDim==2 && srcSpaceDim==2)
563 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
564 INTERP_KERNEL::Interpolation2D interpolation(*this);
565 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
567 else if(srcMeshDim==3 && srcSpaceDim==3)
569 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
570 INTERP_KERNEL::Interpolation3D interpolation(*this);
571 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
573 else if(srcMeshDim==2 && srcSpaceDim==3)
575 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
576 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
577 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
580 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
582 else if(srcMeshDim==-1)
584 if(trgMeshDim==2 && trgSpaceDim==2)
586 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
587 INTERP_KERNEL::Interpolation2D interpolation(*this);
588 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
590 else if(trgMeshDim==3 && trgSpaceDim==3)
592 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
593 INTERP_KERNEL::Interpolation3D interpolation(*this);
594 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
596 else if(trgMeshDim==2 && trgSpaceDim==3)
598 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
599 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
600 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
603 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
606 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
607 _deno_multiply.clear();
608 _deno_multiply.resize(_matrix.size());
609 _deno_reverse_multiply.clear();
610 _deno_reverse_multiply.resize(nbCols);
615 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
617 std::string srcMeth,trgMeth;
618 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
619 const MEDCouplingMappedExtrudedMesh *src_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_src_ft->getMesh());
620 const MEDCouplingMappedExtrudedMesh *target_mesh=static_cast<const MEDCouplingMappedExtrudedMesh *>(_target_ft->getMesh());
622 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
623 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
624 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
625 INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
626 std::vector<std::map<int,double> > matrix2D;
627 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
628 MEDCouplingUMesh *s1D,*t1D;
630 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
631 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
632 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
633 std::vector<std::map<int,double> > matrix1D;
634 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
635 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
638 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
639 target_mesh->getMesh3DIds()->getConstPointer());
641 _deno_multiply.clear();
642 _deno_multiply.resize(_matrix.size());
643 _deno_reverse_multiply.clear();
644 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
649 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
651 std::string srcMeth,trgMeth;
652 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
653 if(methodCpp!="P0P0")
654 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
655 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
656 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
657 const int srcMeshDim=src_mesh->getMeshDimension();
658 const int srcSpceDim=src_mesh->getSpaceDimension();
659 const int trgMeshDim=target_mesh->getMeshDimension();
660 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
661 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 !");
662 std::vector<std::map<int,double> > res;
667 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
668 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
669 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
670 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
675 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
676 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
677 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
678 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
683 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
684 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
685 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
686 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
690 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
692 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
693 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
695 _deno_multiply.clear();
696 _deno_multiply.resize(_matrix.size());
697 _deno_reverse_multiply.clear();
698 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
703 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
705 std::string srcMeth,trgMeth;
706 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
707 if(methodCpp!="P0P0")
708 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
709 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
710 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
711 const int srcMeshDim=src_mesh->getMeshDimension();
712 const int trgMeshDim=target_mesh->getMeshDimension();
713 const int trgSpceDim=target_mesh->getSpaceDimension();
714 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
715 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 !");
720 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
721 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
722 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
723 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
728 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
729 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
730 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
731 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
736 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
737 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
738 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
739 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
743 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
745 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
747 _deno_multiply.clear();
748 _deno_multiply.resize(_matrix.size());
749 _deno_reverse_multiply.clear();
750 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
755 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
757 std::string srcMeth,trgMeth;
758 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
759 if(methodCpp!="P0P0")
760 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
761 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
762 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
763 const int srcMeshDim=src_mesh->getMeshDimension();
764 const int trgMeshDim=target_mesh->getMeshDimension();
765 if(trgMeshDim!=srcMeshDim)
766 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
771 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
772 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
773 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
774 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
779 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
780 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
781 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
782 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
787 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
788 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
789 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
790 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
794 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
796 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
798 _deno_multiply.clear();
799 _deno_multiply.resize(_matrix.size());
800 _deno_reverse_multiply.clear();
801 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
806 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
808 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
809 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 !");
810 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
811 const double *trgLocPtr=trgLoc->begin();
812 int trgSpaceDim=trgLoc->getNumberOfComponents();
813 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
814 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
816 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
817 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
818 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
819 throw INTERP_KERNEL::Exception(oss.str().c_str());
821 const int *srcOffsetArrPtr=srcOffsetArr->begin();
822 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
823 const double *srcLocPtr=srcLoc->begin();
824 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
825 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
826 _matrix.resize(trgNbOfGaussPts);
827 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
828 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
829 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
830 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
831 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
833 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
834 int srcCellId=elts[eltsIndex[*trgId]];
835 double dist=std::numeric_limits<double>::max();
837 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
839 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
841 for(int i=0;i<trgSpaceDim;i++)
842 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
844 { dist=tmp; srcEntry=srcId; }
846 _matrix[*trgId][srcEntry]=1.;
848 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
850 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
851 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
852 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
853 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
854 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
855 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
857 _deno_multiply.clear();
858 _deno_multiply.resize(_matrix.size());
859 _deno_reverse_multiply.clear();
860 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
866 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
867 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
869 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
871 if(method=="GAUSSGAUSS")
873 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
874 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
875 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
876 throw INTERP_KERNEL::Exception(oss.str().c_str());
880 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
881 * 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
882 * only INTERP_KERNEL method should be applied.
884 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
886 std::string srcm,trgm,method;
887 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
888 switch(_interp_matrix_pol)
890 case IK_ONLY_PREFERED:
894 std::string tmp1,tmp2;
895 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
898 catch(INTERP_KERNEL::Exception& /*e*/)
903 case NOT_IK_ONLY_PREFERED:
907 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
910 catch(INTERP_KERNEL::Exception& /*e*/)
917 case NOT_IK_ONLY_FORCED:
920 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
924 void MEDCouplingRemapper::updateTime() const
928 void MEDCouplingRemapper::checkPrepare() const
930 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
932 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
933 if(!s->getMesh() || !t->getMesh())
934 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
938 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
939 * This method returns 3 informations (2 in ouput parameters and 1 in return).
941 * \param [out] srcMeth the string code of the discretization of source field template
942 * \param [out] trgMeth the string code of the discretization of target field template
943 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
945 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
947 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
949 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
950 if(!s->getMesh() || !t->getMesh())
951 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
952 srcMeth=_src_ft->getDiscretization()->getRepr();
953 trgMeth=_target_ft->getDiscretization()->getRepr();
954 return BuildMethodFrom(srcMeth,trgMeth);
957 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
959 std::string method(meth1); method+=meth2;
963 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
967 if(matrixSuppression)
970 _deno_multiply.clear();
971 _deno_reverse_multiply.clear();
975 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
977 if(!srcField || !targetField)
978 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
979 srcField->checkConsistencyLight();
981 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
982 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
983 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
984 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
985 if(srcField->getNature()!=targetField->getNature())
986 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
987 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
989 std::ostringstream oss;
990 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();
991 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
992 throw INTERP_KERNEL::Exception(oss.str().c_str());
994 DataArrayDouble *array(targetField->getArray());
995 int srcNbOfCompo(srcField->getNumberOfComponents());
998 targetField->checkConsistencyLight();
999 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1000 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1005 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1006 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1007 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1008 targetField->setArray(tmp);
1010 computeDeno(srcField->getNature(),srcField,targetField);
1011 double *resPointer(targetField->getArray()->getPointer());
1012 const double *inputPointer(srcField->getArray()->getConstPointer());
1013 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1016 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1019 return computeDenoFromScratch(nat,srcField,trgField);
1020 else if(nat!=_nature_of_deno)
1021 return computeDenoFromScratch(nat,srcField,trgField);
1022 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1023 return computeDenoFromScratch(nat,srcField,trgField);
1026 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1028 _nature_of_deno=nat;
1029 _time_deno_update=getTimeOfThis();
1030 switch(_nature_of_deno)
1032 case IntensiveMaximum:
1034 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1037 case ExtensiveMaximum:
1039 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1040 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1041 const double *denoPtr=deno->getArray()->getConstPointer();
1042 const double *denoRPtr=denoR->getArray()->getConstPointer();
1043 if(trgField->getMesh()->getMeshDimension()==-1)
1045 double *denoRPtr2=denoR->getArray()->getPointer();
1046 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1048 if(srcField->getMesh()->getMeshDimension()==-1)
1050 double *denoPtr2=deno->getArray()->getPointer();
1051 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1054 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1055 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1057 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1058 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1064 case ExtensiveConservation:
1066 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1069 case IntensiveConservation:
1071 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1072 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1073 const double *denoPtr=deno->getArray()->getConstPointer();
1074 const double *denoRPtr=denoR->getArray()->getConstPointer();
1075 if(trgField->getMesh()->getMeshDimension()==-1)
1077 double *denoRPtr2=denoR->getArray()->getPointer();
1078 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1080 if(srcField->getMesh()->getMeshDimension()==-1)
1082 double *denoPtr2=deno->getArray()->getPointer();
1083 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1086 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1087 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1089 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1090 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1097 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1101 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1104 double *tmp=new double[inputNbOfCompo];
1105 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1107 if((*iter1).empty())
1110 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1114 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1115 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1116 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1118 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1119 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1125 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1127 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1129 double *tmp=new double[inputNbOfCompo];
1130 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1131 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1133 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1135 isReached[(*iter2).first]=true;
1136 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1137 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1142 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1144 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1147 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1149 matOut.resize(nbColsMatIn);
1151 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1152 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1153 matOut[(*iter2).first][id]=(*iter2).second;
1156 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1157 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1159 std::map<int,double> values;
1161 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1164 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1166 sum+=(*iter2).second;
1167 values[(*iter2).first]+=(*iter2).second;
1169 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1170 deno[idx][(*iter2).first]=sum;
1173 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1175 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1176 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1180 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1181 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1183 std::map<int,double> values;
1185 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1188 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1190 sum+=(*iter2).second;
1191 values[(*iter2).first]+=(*iter2).second;
1193 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1194 denoReverse[(*iter2).first][idx]=sum;
1197 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1199 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1200 deno[idx][(*iter2).first]=values[(*iter2).first];
1204 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1205 const std::vector< std::map<int,double> >& m2D,
1206 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1207 const int *corrCellIdTrg)
1209 int nbOf2DCellsTrg=m2D.size();
1210 int nbOf1DCellsTrg=m1D.size();
1211 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1212 _matrix.resize(nbOf3DCellsTrg);
1214 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1216 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1219 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1221 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1223 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1230 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1233 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1235 std::cout << "Target Cell # " << id << " : ";
1236 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1237 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1238 std::cout << std::endl;
1242 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1248 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1250 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1252 return (int)_deno_reverse_multiply.size();
1256 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1257 * If not the behaviour is unpredictable.
1258 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1259 * 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.
1260 * 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
1263 * \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.
1264 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1265 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1267 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1270 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1272 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1274 std::map<int,double>& rowNew=matrixNew[i];
1275 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1277 if(fabs((*it2).second)>maxValAbs)
1278 rowNew[(*it2).first]=(*it2).second;
1289 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1290 * If not the behaviour is unpredictable.
1291 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1292 * 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.
1293 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1294 * 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
1297 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1298 * \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
1299 * that all coefficients are null.
1300 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1302 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1304 double maxVal=getMaxValueInCrudeMatrix();
1307 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1311 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1312 * If not the behaviour is unpredictable.
1313 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1314 * The returned value is positive.
1316 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1319 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1320 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1321 if(fabs((*it2).second)>ret)
1322 ret=fabs((*it2).second);