1 // Copyright (C) 2007-2013 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.
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 "MEDCouplingExtrudedMesh.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 "Interpolation3D2D.txx"
38 #include "InterpolationCU.txx"
39 #include "InterpolationCC.txx"
41 using namespace ParaMEDMEM;
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()
52 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const std::string& method)
54 if(!srcMesh || !targetMesh)
55 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepare : presence of NULL input pointer !");
56 std::string srcMethod,targetMethod;
57 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,srcMethod,targetMethod);
58 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldTemplate> src=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(srcMethod));
59 src->setMesh(srcMesh);
60 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldTemplate> target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod));
61 target->setMesh(targetMesh);
62 return prepareEx(src,target);
65 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target)
68 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL input pointer !");
69 if(!src->getMesh() || !target->getMesh())
70 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareEx : presence of NULL mesh pointer in given field template !");
72 _src_ft=const_cast<MEDCouplingFieldTemplate *>(src); _src_ft->incrRef();
73 _target_ft=const_cast<MEDCouplingFieldTemplate *>(target); _target_ft->incrRef();
74 if(isInterpKernelOnlyOrNotOnly())
75 return prepareInterpKernelOnly();
77 return prepareNotInterpKernelOnly();
80 int MEDCouplingRemapper::prepareInterpKernelOnly()
82 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
83 switch(meshInterpType)
93 case 85://Unstructured-Unstructured
94 return prepareInterpKernelOnlyUU();
97 case 87://Unstructured-Cartesian
98 return prepareInterpKernelOnlyUC();
101 case 117://Cartesian-Unstructured
102 return prepareInterpKernelOnlyCU();
103 case 119://Cartesian-Cartesian
104 return prepareInterpKernelOnlyCC();
105 case 136://Extruded-Extruded
106 return prepareInterpKernelOnlyEE();
108 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
112 int MEDCouplingRemapper::prepareNotInterpKernelOnly()
114 std::string srcm,trgm,method;
115 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
116 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
119 return prepareNotInterpKernelOnlyGaussGauss();
122 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
123 throw INTERP_KERNEL::Exception(oss.str().c_str());
129 * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
130 * If meshes of \b srcField and \b targetField do not match exactly those given into \ref ParaMEDMEM::MEDCouplingRemapper::prepare "prepare method" an exception will be thrown.
132 * \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 ParaMEDMEM::MEDCouplingRemapper::prepare.
133 * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
135 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue)
137 transferUnderground(srcField,targetField,true,dftValue);
141 * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
142 * 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
144 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
146 * \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 ParaMEDMEM::MEDCouplingRemapper::prepare.
147 * \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.
149 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField)
151 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
154 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue)
157 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
158 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
159 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
160 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
161 if(srcField->getNature()!=targetField->getNature())
162 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
163 DataArrayDouble *array=srcField->getArray();
164 int trgNbOfCompo=targetField->getNumberOfComponents();
167 if(trgNbOfCompo!=srcField->getNumberOfComponents())
168 throw INTERP_KERNEL::Exception("Number of components mismatch !");
172 array=DataArrayDouble::New();
173 array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
174 srcField->setArray(array);
177 computeDeno(srcField->getNature(),srcField,targetField);
178 double *resPointer=array->getPointer();
179 const double *inputPointer=targetField->getArray()->getConstPointer();
180 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
183 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue)
186 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
187 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
188 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
189 ret->setNature(srcField->getNature());
190 transfer(srcField,ret,dftValue);
191 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
195 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue)
198 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
199 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
200 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
201 ret->setNature(targetField->getNature());
202 reverseTransfer(ret,targetField,dftValue);
203 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
208 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
209 * is here only for automatic CORBA generators.
211 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
213 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
217 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
218 * is here only for automatic CORBA generators.
220 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
222 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
226 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
227 * is here only for automatic CORBA generators.
229 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
231 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
235 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or prefered.
236 * If interpolation matrix policy is :
238 * - 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
239 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
240 * If not, the \b not only INTERP_KERNEL method will be attempt.
242 * - 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
243 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
244 * If not, the INTERP_KERNEL only method will be attempt.
246 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
248 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
250 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
252 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
254 return _interp_matrix_pol;
258 * 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
259 * 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
260 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
262 * If interpolation matrix policy is :
264 * - 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
265 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
266 * If not, the \b not only INTERP_KERNEL method will be attempt.
268 * - 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
269 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
270 * If not, the INTERP_KERNEL only method will be attempt.
272 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
274 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
276 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c ParaMEDMEM::InterpolationMatrixPolicy
277 * for automatic generation of CORBA component.
279 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
281 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol)
283 switch(newInterpMatPol)
286 _interp_matrix_pol=IK_ONLY_PREFERED;
289 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
292 _interp_matrix_pol=IK_ONLY_FORCED;
295 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
298 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 !");
302 int MEDCouplingRemapper::prepareInterpKernelOnlyUU()
304 const MEDCouplingPointSet *src_mesh=static_cast<const MEDCouplingPointSet *>(_src_ft->getMesh());
305 const MEDCouplingPointSet *target_mesh=static_cast<const MEDCouplingPointSet *>(_target_ft->getMesh());
306 std::string srcMeth,trgMeth;
307 std::string method(checkAndGiveInterpolationMethodStr(srcMeth,trgMeth));
308 const int srcMeshDim=src_mesh->getMeshDimension();
311 srcSpaceDim=src_mesh->getSpaceDimension();
312 const int trgMeshDim=target_mesh->getMeshDimension();
315 trgSpaceDim=target_mesh->getSpaceDimension();
316 if(trgSpaceDim!=srcSpaceDim)
317 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
318 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
320 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
322 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
323 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
324 INTERP_KERNEL::Interpolation1D interpolation(*this);
325 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
327 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
329 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
330 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
331 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
332 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
334 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
336 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
337 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
338 INTERP_KERNEL::Interpolation2D interpolation(*this);
339 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
341 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
343 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
344 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
345 INTERP_KERNEL::Interpolation3D interpolation(*this);
346 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
348 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
350 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
351 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
352 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
353 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
355 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
357 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
358 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
359 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
360 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
361 INTERP_KERNEL::Interpolation3D interpolation(*this);
362 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
364 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
366 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
367 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
368 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
369 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
370 INTERP_KERNEL::Interpolation3D interpolation(*this);
371 std::vector<std::map<int,double> > matrixTmp;
372 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
373 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
374 ReverseMatrix(matrixTmp,nbCols,_matrix);
375 nbCols=matrixTmp.size();
377 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
379 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
381 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
382 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
383 INTERP_KERNEL::Interpolation2D interpolation(*this);
384 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
388 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
389 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
390 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
391 std::vector<std::map<int,double> > matrixTmp;
392 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
393 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
394 ReverseMatrix(matrixTmp,nbCols,_matrix);
395 nbCols=matrixTmp.size();
396 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
397 if(!duplicateFaces.empty())
399 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
400 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
402 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
403 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
409 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
411 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
413 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
414 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
415 INTERP_KERNEL::Interpolation2D interpolation(*this);
416 std::vector<std::map<int,double> > matrixTmp;
417 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
418 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
419 ReverseMatrix(matrixTmp,nbCols,_matrix);
420 nbCols=matrixTmp.size();
424 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
425 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
426 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
427 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
428 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
429 if(!duplicateFaces.empty())
431 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
432 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
434 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
435 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
441 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
443 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
444 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
445 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
446 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
447 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
448 if(!duplicateFaces.empty())
450 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
451 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
453 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
454 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
459 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
461 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
462 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
463 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
464 std::vector<std::map<int,double> > matrixTmp;
465 std::string revMethod(BuildMethodFrom(trgMeth,srcMeth));
466 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,revMethod);
467 ReverseMatrix(matrixTmp,nbCols,_matrix);
468 nbCols=matrixTmp.size();
469 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
470 if(!duplicateFaces.empty())
472 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
473 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
475 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
476 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
481 else if(trgMeshDim==-1)
483 if(srcMeshDim==2 && srcSpaceDim==2)
485 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
486 INTERP_KERNEL::Interpolation2D interpolation(*this);
487 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
489 else if(srcMeshDim==3 && srcSpaceDim==3)
491 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
492 INTERP_KERNEL::Interpolation3D interpolation(*this);
493 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
495 else if(srcMeshDim==2 && srcSpaceDim==3)
497 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
498 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
499 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth);
502 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
504 else if(srcMeshDim==-1)
506 if(trgMeshDim==2 && trgSpaceDim==2)
508 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
509 INTERP_KERNEL::Interpolation2D interpolation(*this);
510 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
512 else if(trgMeshDim==3 && trgSpaceDim==3)
514 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
515 INTERP_KERNEL::Interpolation3D interpolation(*this);
516 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
518 else if(trgMeshDim==2 && trgSpaceDim==3)
520 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
521 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
522 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth);
525 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
528 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
529 _deno_multiply.clear();
530 _deno_multiply.resize(_matrix.size());
531 _deno_reverse_multiply.clear();
532 _deno_reverse_multiply.resize(nbCols);
537 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
539 std::string srcMeth,trgMeth;
540 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
541 const MEDCouplingExtrudedMesh *src_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_src_ft->getMesh());
542 const MEDCouplingExtrudedMesh *target_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_target_ft->getMesh());
544 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
545 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
546 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
547 INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
548 std::vector<std::map<int,double> > matrix2D;
549 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC);
550 MEDCouplingUMesh *s1D,*t1D;
552 MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
553 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
554 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
555 std::vector<std::map<int,double> > matrix1D;
556 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
557 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC);
560 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
561 target_mesh->getMesh3DIds()->getConstPointer());
563 _deno_multiply.clear();
564 _deno_multiply.resize(_matrix.size());
565 _deno_reverse_multiply.clear();
566 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
571 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
573 std::string srcMeth,trgMeth;
574 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
575 if(methodCpp!="P0P0")
576 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
577 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
578 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
579 const int srcMeshDim=src_mesh->getMeshDimension();
580 const int srcSpceDim=src_mesh->getSpaceDimension();
581 const int trgMeshDim=target_mesh->getMeshDimension();
582 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
583 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 !");
584 std::vector<std::map<int,double> > res;
589 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
590 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
591 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
592 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
597 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
598 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
599 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
600 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
605 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
606 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
607 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
608 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
612 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
614 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
615 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
617 _deno_multiply.clear();
618 _deno_multiply.resize(_matrix.size());
619 _deno_reverse_multiply.clear();
620 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
625 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
627 std::string srcMeth,trgMeth;
628 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
629 if(methodCpp!="P0P0")
630 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
631 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
632 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
633 const int srcMeshDim=src_mesh->getMeshDimension();
634 const int trgMeshDim=target_mesh->getMeshDimension();
635 const int trgSpceDim=target_mesh->getSpaceDimension();
636 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
637 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 !");
642 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
643 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
644 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
645 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
650 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
651 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
652 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
653 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
658 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
659 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
660 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
661 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
665 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
667 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
669 _deno_multiply.clear();
670 _deno_multiply.resize(_matrix.size());
671 _deno_reverse_multiply.clear();
672 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
677 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
679 std::string srcMeth,trgMeth;
680 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
681 if(methodCpp!="P0P0")
682 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
683 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
684 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
685 const int srcMeshDim=src_mesh->getMeshDimension();
686 const int trgMeshDim=target_mesh->getMeshDimension();
687 if(trgMeshDim!=srcMeshDim)
688 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
693 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
694 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
695 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
696 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
701 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
702 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
703 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
704 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
709 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
710 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
711 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
712 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
716 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
718 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
720 _deno_multiply.clear();
721 _deno_multiply.resize(_matrix.size());
722 _deno_reverse_multiply.clear();
723 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
728 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
730 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
731 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 !");
732 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
733 const double *trgLocPtr=trgLoc->begin();
734 int trgSpaceDim=trgLoc->getNumberOfComponents();
735 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
736 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
738 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
739 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
740 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
741 throw INTERP_KERNEL::Exception(oss.str().c_str());
743 const int *srcOffsetArrPtr=srcOffsetArr->begin();
744 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
745 const double *srcLocPtr=srcLoc->begin();
746 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsArr,eltsIndexArr;
747 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
748 _matrix.resize(trgNbOfGaussPts);
749 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
750 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
751 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
752 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->getIdsNotEqual(0);
753 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
755 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
756 int srcCellId=elts[eltsIndex[*trgId]];
757 double dist=std::numeric_limits<double>::max();
759 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
761 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
763 for(int i=0;i<trgSpaceDim;i++)
764 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
766 { dist=tmp; srcEntry=srcId; }
768 _matrix[*trgId][srcEntry]=1.;
770 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
772 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->getIdsEqual(0);
773 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
774 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
775 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
776 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
777 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
779 _deno_multiply.clear();
780 _deno_multiply.resize(_matrix.size());
781 _deno_reverse_multiply.clear();
782 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
788 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
789 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
791 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
793 if(method=="GAUSSGAUSS")
795 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
796 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
797 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
798 throw INTERP_KERNEL::Exception(oss.str().c_str());
802 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
803 * to IK_ONLY_PREFERED = 0 ) , which method will be applied. If \c true is returned the INTERP_KERNEL only method should be applied to \c false the \b not
804 * only INTERP_KERNEL method should be applied.
806 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
808 std::string srcm,trgm,method;
809 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
810 switch(_interp_matrix_pol)
812 case IK_ONLY_PREFERED:
816 std::string tmp1,tmp2;
817 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method,tmp1,tmp2);
820 catch(INTERP_KERNEL::Exception& /*e*/)
825 case NOT_IK_ONLY_PREFERED:
829 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
832 catch(INTERP_KERNEL::Exception& /*e*/)
839 case NOT_IK_ONLY_FORCED:
842 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
846 void MEDCouplingRemapper::updateTime() const
850 void MEDCouplingRemapper::checkPrepare() const
852 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
854 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
855 if(!s->getMesh() || !t->getMesh())
856 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
860 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
861 * This method returns 3 informations (2 in ouput parameters and 1 in return).
863 * \param [out] srcMeth the string code of the discretization of source field template
864 * \param [out] trgMeth the string code of the discretization of target field template
865 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
867 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
869 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
871 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
872 if(!s->getMesh() || !t->getMesh())
873 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
874 srcMeth=_src_ft->getDiscretization()->getRepr();
875 trgMeth=_target_ft->getDiscretization()->getRepr();
876 return BuildMethodFrom(srcMeth,trgMeth);
879 std::string MEDCouplingRemapper::BuildMethodFrom(const std::string& meth1, const std::string& meth2)
881 std::string method(meth1); method+=meth2;
885 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
889 if(matrixSuppression)
892 _deno_multiply.clear();
893 _deno_reverse_multiply.clear();
897 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
900 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
901 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
902 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
903 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
904 if(srcField->getNature()!=targetField->getNature())
905 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
906 DataArrayDouble *array=targetField->getArray();
907 int srcNbOfCompo=srcField->getNumberOfComponents();
910 if(srcNbOfCompo!=targetField->getNumberOfComponents())
911 throw INTERP_KERNEL::Exception("Number of components mismatch !");
916 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
917 array=DataArrayDouble::New();
918 array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
919 targetField->setArray(array);
922 computeDeno(srcField->getNature(),srcField,targetField);
923 double *resPointer=array->getPointer();
924 const double *inputPointer=srcField->getArray()->getConstPointer();
925 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
928 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
931 return computeDenoFromScratch(nat,srcField,trgField);
932 else if(nat!=_nature_of_deno)
933 return computeDenoFromScratch(nat,srcField,trgField);
934 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
935 return computeDenoFromScratch(nat,srcField,trgField);
938 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
941 _time_deno_update=getTimeOfThis();
942 switch(_nature_of_deno)
944 case ConservativeVolumic:
946 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
951 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
952 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
953 const double *denoPtr=deno->getArray()->getConstPointer();
954 const double *denoRPtr=denoR->getArray()->getConstPointer();
955 if(trgField->getMesh()->getMeshDimension()==-1)
957 double *denoRPtr2=denoR->getArray()->getPointer();
958 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
960 if(srcField->getMesh()->getMeshDimension()==-1)
962 double *denoPtr2=deno->getArray()->getPointer();
963 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
966 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
967 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
969 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
970 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
976 case IntegralGlobConstraint:
978 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
983 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
984 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
985 const double *denoPtr=deno->getArray()->getConstPointer();
986 const double *denoRPtr=denoR->getArray()->getConstPointer();
987 if(trgField->getMesh()->getMeshDimension()==-1)
989 double *denoRPtr2=denoR->getArray()->getPointer();
990 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
992 if(srcField->getMesh()->getMeshDimension()==-1)
994 double *denoPtr2=deno->getArray()->getPointer();
995 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
998 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
999 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1001 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
1002 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1009 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1013 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1016 double *tmp=new double[inputNbOfCompo];
1017 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1019 if((*iter1).empty())
1022 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1026 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1027 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1028 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1030 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1031 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1037 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1039 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1041 double *tmp=new double[inputNbOfCompo];
1042 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1043 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1045 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1047 isReached[(*iter2).first]=true;
1048 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1049 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1054 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1056 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1059 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1061 matOut.resize(nbColsMatIn);
1063 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1064 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1065 matOut[(*iter2).first][id]=(*iter2).second;
1068 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1069 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1071 std::map<int,double> values;
1073 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1076 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1078 sum+=(*iter2).second;
1079 values[(*iter2).first]+=(*iter2).second;
1081 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1082 deno[idx][(*iter2).first]=sum;
1085 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1087 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1088 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1092 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1093 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1095 std::map<int,double> values;
1097 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1100 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1102 sum+=(*iter2).second;
1103 values[(*iter2).first]+=(*iter2).second;
1105 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1106 denoReverse[(*iter2).first][idx]=sum;
1109 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1111 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1112 deno[idx][(*iter2).first]=values[(*iter2).first];
1116 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1117 const std::vector< std::map<int,double> >& m2D,
1118 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1119 const int *corrCellIdTrg)
1121 int nbOf2DCellsTrg=m2D.size();
1122 int nbOf1DCellsTrg=m1D.size();
1123 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1124 _matrix.resize(nbOf3DCellsTrg);
1126 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1128 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1131 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1133 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1135 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1142 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1145 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1147 std::cout << "Target Cell # " << id << " : ";
1148 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1149 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1150 std::cout << std::endl;
1154 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1160 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1162 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1164 return (int)_deno_reverse_multiply.size();
1168 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1169 * If not the behaviour is unpredictable.
1170 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1171 * 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.
1172 * 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
1175 * \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.
1176 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1177 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1179 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1182 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1184 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1186 std::map<int,double>& rowNew=matrixNew[i];
1187 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1189 if(fabs((*it2).second)>maxValAbs)
1190 rowNew[(*it2).first]=(*it2).second;
1201 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1202 * If not the behaviour is unpredictable.
1203 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1204 * 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.
1205 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1206 * 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
1209 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1210 * \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
1211 * that all coefficients are null.
1212 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1214 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1216 double maxVal=getMaxValueInCrudeMatrix();
1219 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1223 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1224 * If not the behaviour is unpredictable.
1225 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1226 * The returned value is positive.
1228 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1231 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1232 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1233 if(fabs((*it2).second)>ret)
1234 ret=fabs((*it2).second);