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 MCAuto<MEDCouplingUMesh> src2D(src_mesh->getMesh2D()->clone(false)); src2D->changeSpaceDimension(2,0.);
624 MCAuto<MEDCouplingUMesh> trg2D(target_mesh->getMesh2D()->clone(false)); trg2D->changeSpaceDimension(2,0.);
625 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src2D);
626 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(trg2D);
627 INTERP_KERNEL::Interpolation2D interpolation2D(*this);
628 std::vector<std::map<int,double> > matrix2D;
629 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
630 MEDCouplingUMesh *s1D,*t1D;
632 MEDCouplingMappedExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
633 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
634 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
635 std::vector<std::map<int,double> > matrix1D;
636 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
637 if(interpolation1D.getIntersectionType()==INTERP_KERNEL::Geometric2D)// For intersection type of 1D, Geometric2D do not deal with it ! -> make interpolation1D not inherite from this
638 interpolation1D.setIntersectionType(INTERP_KERNEL::Triangulation);//
639 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
642 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
643 target_mesh->getMesh3DIds()->getConstPointer());
645 _deno_multiply.clear();
646 _deno_multiply.resize(_matrix.size());
647 _deno_reverse_multiply.clear();
648 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
653 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
655 std::string srcMeth,trgMeth;
656 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
657 if(methodCpp!="P0P0")
658 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
659 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
660 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
661 const int srcMeshDim=src_mesh->getMeshDimension();
662 const int srcSpceDim=src_mesh->getSpaceDimension();
663 const int trgMeshDim=target_mesh->getMeshDimension();
664 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
665 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 !");
666 std::vector<std::map<int,double> > res;
671 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
672 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
673 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
674 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
679 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
680 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
681 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
682 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
687 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
688 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
689 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
690 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
694 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
696 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
697 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
699 _deno_multiply.clear();
700 _deno_multiply.resize(_matrix.size());
701 _deno_reverse_multiply.clear();
702 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
707 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
709 std::string srcMeth,trgMeth;
710 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
711 if(methodCpp!="P0P0")
712 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
713 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
714 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
715 const int srcMeshDim=src_mesh->getMeshDimension();
716 const int trgMeshDim=target_mesh->getMeshDimension();
717 const int trgSpceDim=target_mesh->getSpaceDimension();
718 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
719 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 !");
724 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
725 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
726 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
727 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
732 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
733 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
734 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
735 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
740 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
741 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
742 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
743 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
747 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
749 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
751 _deno_multiply.clear();
752 _deno_multiply.resize(_matrix.size());
753 _deno_reverse_multiply.clear();
754 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
759 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
761 std::string srcMeth,trgMeth;
762 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
763 if(methodCpp!="P0P0")
764 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
765 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
766 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
767 const int srcMeshDim=src_mesh->getMeshDimension();
768 const int trgMeshDim=target_mesh->getMeshDimension();
769 if(trgMeshDim!=srcMeshDim)
770 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
775 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
776 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
777 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
778 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
783 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
784 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
785 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
786 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
791 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
792 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
793 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
794 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
798 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
800 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
802 _deno_multiply.clear();
803 _deno_multiply.resize(_matrix.size());
804 _deno_reverse_multiply.clear();
805 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
810 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
812 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
813 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 !");
814 MCAuto<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
815 const double *trgLocPtr=trgLoc->begin();
816 int trgSpaceDim=trgLoc->getNumberOfComponents();
817 MCAuto<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
818 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
820 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
821 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
822 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
823 throw INTERP_KERNEL::Exception(oss.str().c_str());
825 const int *srcOffsetArrPtr=srcOffsetArr->begin();
826 MCAuto<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
827 const double *srcLocPtr=srcLoc->begin();
828 MCAuto<DataArrayInt> eltsArr,eltsIndexArr;
829 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
830 _matrix.resize(trgNbOfGaussPts);
831 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
832 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
833 MCAuto<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
834 MCAuto<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->findIdsNotEqual(0);
835 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
837 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
838 int srcCellId=elts[eltsIndex[*trgId]];
839 double dist=std::numeric_limits<double>::max();
841 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
843 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
845 for(int i=0;i<trgSpaceDim;i++)
846 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
848 { dist=tmp; srcEntry=srcId; }
850 _matrix[*trgId][srcEntry]=1.;
852 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
854 MCAuto<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->findIdsEqual(0);
855 MCAuto<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
856 MCAuto<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
857 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
858 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
859 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
861 _deno_multiply.clear();
862 _deno_multiply.resize(_matrix.size());
863 _deno_reverse_multiply.clear();
864 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
870 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
871 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
873 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
875 if(method=="GAUSSGAUSS")
877 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
878 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
879 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
880 throw INTERP_KERNEL::Exception(oss.str().c_str());
884 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
885 * 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
886 * only INTERP_KERNEL method should be applied.
888 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
890 std::string srcm,trgm,method;
891 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
892 switch(_interp_matrix_pol)
894 case IK_ONLY_PREFERED:
898 std::string tmp1,tmp2;
899 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
902 catch(INTERP_KERNEL::Exception& /*e*/)
907 case NOT_IK_ONLY_PREFERED:
911 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
914 catch(INTERP_KERNEL::Exception& /*e*/)
921 case NOT_IK_ONLY_FORCED:
924 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
928 void MEDCouplingRemapper::updateTime() const
932 void MEDCouplingRemapper::checkPrepare() const
934 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
936 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
937 if(!s->getMesh() || !t->getMesh())
938 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
942 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
943 * This method returns 3 informations (2 in ouput parameters and 1 in return).
945 * \param [out] srcMeth the string code of the discretization of source field template
946 * \param [out] trgMeth the string code of the discretization of target field template
947 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
949 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
951 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
953 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
954 if(!s->getMesh() || !t->getMesh())
955 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
956 srcMeth=_src_ft->getDiscretization()->getRepr();
957 trgMeth=_target_ft->getDiscretization()->getRepr();
958 return BuildMethodFrom(srcMeth,trgMeth);
961 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
963 std::string method(meth1); method+=meth2;
967 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
971 if(matrixSuppression)
974 _deno_multiply.clear();
975 _deno_reverse_multiply.clear();
979 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
981 if(!srcField || !targetField)
982 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::transferUnderground : srcField or targetField is NULL !");
983 srcField->checkConsistencyLight();
985 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
986 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
987 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
988 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
989 if(srcField->getNature()!=targetField->getNature())
990 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
991 if(srcField->getNumberOfTuplesExpected()!=_src_ft->getNumberOfTuplesExpected())
993 std::ostringstream oss;
994 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();
995 oss << " ! It appears that the source support is not the same between the prepare and the transfer !";
996 throw INTERP_KERNEL::Exception(oss.str().c_str());
998 DataArrayDouble *array(targetField->getArray());
999 int srcNbOfCompo(srcField->getNumberOfComponents());
1002 targetField->checkConsistencyLight();
1003 if(srcNbOfCompo!=targetField->getNumberOfComponents())
1004 throw INTERP_KERNEL::Exception("Number of components mismatch !");
1009 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
1010 MCAuto<DataArrayDouble> tmp(DataArrayDouble::New());
1011 tmp->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
1012 targetField->setArray(tmp);
1014 computeDeno(srcField->getNature(),srcField,targetField);
1015 double *resPointer(targetField->getArray()->getPointer());
1016 const double *inputPointer(srcField->getArray()->getConstPointer());
1017 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
1020 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1023 return computeDenoFromScratch(nat,srcField,trgField);
1024 else if(nat!=_nature_of_deno)
1025 return computeDenoFromScratch(nat,srcField,trgField);
1026 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
1027 return computeDenoFromScratch(nat,srcField,trgField);
1030 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
1032 _nature_of_deno=nat;
1033 _time_deno_update=getTimeOfThis();
1034 switch(_nature_of_deno)
1036 case IntensiveMaximum:
1038 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1041 case ExtensiveMaximum:
1043 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1044 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1045 const double *denoPtr=deno->getArray()->getConstPointer();
1046 const double *denoRPtr=denoR->getArray()->getConstPointer();
1047 if(trgField->getMesh()->getMeshDimension()==-1)
1049 double *denoRPtr2=denoR->getArray()->getPointer();
1050 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1052 if(srcField->getMesh()->getMeshDimension()==-1)
1054 double *denoPtr2=deno->getArray()->getPointer();
1055 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1058 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1059 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1061 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
1062 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
1068 case ExtensiveConservation:
1070 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
1073 case IntensiveConservation:
1075 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
1076 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
1077 const double *denoPtr=deno->getArray()->getConstPointer();
1078 const double *denoRPtr=denoR->getArray()->getConstPointer();
1079 if(trgField->getMesh()->getMeshDimension()==-1)
1081 double *denoRPtr2=denoR->getArray()->getPointer();
1082 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
1084 if(srcField->getMesh()->getMeshDimension()==-1)
1086 double *denoPtr2=deno->getArray()->getPointer();
1087 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
1090 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1091 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1093 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1094 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1101 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1105 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1108 double *tmp=new double[inputNbOfCompo];
1109 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1111 if((*iter1).empty())
1114 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1118 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1119 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1120 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1122 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1123 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1129 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1131 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1133 double *tmp=new double[inputNbOfCompo];
1134 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1135 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1137 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1139 isReached[(*iter2).first]=true;
1140 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1141 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1146 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1148 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1151 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1153 matOut.resize(nbColsMatIn);
1155 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1156 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1157 matOut[(*iter2).first][id]=(*iter2).second;
1160 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1161 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1163 std::map<int,double> values;
1165 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1168 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1170 sum+=(*iter2).second;
1171 values[(*iter2).first]+=(*iter2).second;
1173 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1174 deno[idx][(*iter2).first]=sum;
1177 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++)
1180 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1184 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1185 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1187 std::map<int,double> values;
1189 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1192 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1194 sum+=(*iter2).second;
1195 values[(*iter2).first]+=(*iter2).second;
1197 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1198 denoReverse[(*iter2).first][idx]=sum;
1201 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++)
1204 deno[idx][(*iter2).first]=values[(*iter2).first];
1208 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1209 const std::vector< std::map<int,double> >& m2D,
1210 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1211 const int *corrCellIdTrg)
1213 int nbOf2DCellsTrg=m2D.size();
1214 int nbOf1DCellsTrg=m1D.size();
1215 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1216 _matrix.resize(nbOf3DCellsTrg);
1218 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1220 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1223 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1225 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1227 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1234 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1237 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1239 std::cout << "Target Cell # " << id << " : ";
1240 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1241 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1242 std::cout << std::endl;
1246 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1252 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1254 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1256 return (int)_deno_reverse_multiply.size();
1260 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1261 * If not the behaviour is unpredictable.
1262 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1263 * 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.
1264 * 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
1267 * \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.
1268 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1269 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1271 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1274 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1276 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1278 std::map<int,double>& rowNew=matrixNew[i];
1279 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1281 if(fabs((*it2).second)>maxValAbs)
1282 rowNew[(*it2).first]=(*it2).second;
1293 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1294 * If not the behaviour is unpredictable.
1295 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1296 * 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.
1297 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1298 * 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
1301 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1302 * \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
1303 * that all coefficients are null.
1304 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1306 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1308 double maxVal=getMaxValueInCrudeMatrix();
1311 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1315 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1316 * If not the behaviour is unpredictable.
1317 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1318 * The returned value is positive.
1320 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1323 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1324 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1325 if(fabs((*it2).second)>ret)
1326 ret=fabs((*it2).second);