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 char *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.c_str()));
59 src->setMesh(srcMesh);
60 MEDCouplingAutoRefCountObjectPtr<MEDCouplingFieldTemplate> target=MEDCouplingFieldTemplate::New(MEDCouplingFieldDiscretization::GetTypeOfFieldFromStringRepr(targetMethod.c_str()));
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.c_str());
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.c_str());
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.c_str());
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.c_str());
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.c_str());
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.c_str());
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 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
373 ReverseMatrix(matrixTmp,nbCols,_matrix);
374 nbCols=matrixTmp.size();
376 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
378 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
380 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
381 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
382 INTERP_KERNEL::Interpolation2D interpolation(*this);
383 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
387 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
388 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
389 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
390 std::vector<std::map<int,double> > matrixTmp;
391 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
392 ReverseMatrix(matrixTmp,nbCols,_matrix);
393 nbCols=matrixTmp.size();
394 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
395 if(!duplicateFaces.empty())
397 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
398 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
400 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
401 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
407 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
409 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
411 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
412 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
413 INTERP_KERNEL::Interpolation2D interpolation(*this);
414 std::vector<std::map<int,double> > matrixTmp;
415 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
416 ReverseMatrix(matrixTmp,nbCols,_matrix);
417 nbCols=matrixTmp.size();
421 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
422 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
423 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
424 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
425 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
426 if(!duplicateFaces.empty())
428 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
429 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
431 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
432 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
438 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
440 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
441 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
442 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
443 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
444 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
445 if(!duplicateFaces.empty())
447 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
448 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
450 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
451 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
456 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
458 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
459 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
460 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
461 std::vector<std::map<int,double> > matrixTmp;
462 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
463 ReverseMatrix(matrixTmp,nbCols,_matrix);
464 nbCols=matrixTmp.size();
465 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
466 if(!duplicateFaces.empty())
468 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
469 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
471 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
472 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
477 else if(trgMeshDim==-1)
479 if(srcMeshDim==2 && srcSpaceDim==2)
481 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
482 INTERP_KERNEL::Interpolation2D interpolation(*this);
483 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
485 else if(srcMeshDim==3 && srcSpaceDim==3)
487 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
488 INTERP_KERNEL::Interpolation3D interpolation(*this);
489 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
491 else if(srcMeshDim==2 && srcSpaceDim==3)
493 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
494 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
495 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
498 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
500 else if(srcMeshDim==-1)
502 if(trgMeshDim==2 && trgSpaceDim==2)
504 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
505 INTERP_KERNEL::Interpolation2D interpolation(*this);
506 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
508 else if(trgMeshDim==3 && trgSpaceDim==3)
510 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
511 INTERP_KERNEL::Interpolation3D interpolation(*this);
512 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
514 else if(trgMeshDim==2 && trgSpaceDim==3)
516 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
517 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
518 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
521 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
524 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
525 _deno_multiply.clear();
526 _deno_multiply.resize(_matrix.size());
527 _deno_reverse_multiply.clear();
528 _deno_reverse_multiply.resize(nbCols);
533 int MEDCouplingRemapper::prepareInterpKernelOnlyEE()
535 std::string srcMeth,trgMeth;
536 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
537 const MEDCouplingExtrudedMesh *src_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_src_ft->getMesh());
538 const MEDCouplingExtrudedMesh *target_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_target_ft->getMesh());
540 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
541 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
542 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
543 INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
544 std::vector<std::map<int,double> > matrix2D;
545 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC.c_str());
546 MEDCouplingUMesh *s1D,*t1D;
548 MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
549 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
550 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
551 std::vector<std::map<int,double> > matrix1D;
552 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
553 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC.c_str());
556 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
557 target_mesh->getMesh3DIds()->getConstPointer());
559 _deno_multiply.clear();
560 _deno_multiply.resize(_matrix.size());
561 _deno_reverse_multiply.clear();
562 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
567 int MEDCouplingRemapper::prepareInterpKernelOnlyUC()
569 std::string srcMeth,trgMeth;
570 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
571 if(methodCpp!="P0P0")
572 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
573 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
574 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
575 const int srcMeshDim=src_mesh->getMeshDimension();
576 const int srcSpceDim=src_mesh->getSpaceDimension();
577 const int trgMeshDim=target_mesh->getMeshDimension();
578 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
579 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 !");
580 std::vector<std::map<int,double> > res;
585 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
586 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
587 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
588 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
593 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
594 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
595 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
596 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
601 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
602 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
603 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
604 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
608 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
610 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
611 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
613 _deno_multiply.clear();
614 _deno_multiply.resize(_matrix.size());
615 _deno_reverse_multiply.clear();
616 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
621 int MEDCouplingRemapper::prepareInterpKernelOnlyCU()
623 std::string srcMeth,trgMeth;
624 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
625 if(methodCpp!="P0P0")
626 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
627 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
628 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
629 const int srcMeshDim=src_mesh->getMeshDimension();
630 const int trgMeshDim=target_mesh->getMeshDimension();
631 const int trgSpceDim=target_mesh->getSpaceDimension();
632 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
633 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 !");
638 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
639 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
640 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
641 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
646 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
647 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
648 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
649 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
654 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
655 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
656 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
657 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
661 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
663 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
665 _deno_multiply.clear();
666 _deno_multiply.resize(_matrix.size());
667 _deno_reverse_multiply.clear();
668 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
673 int MEDCouplingRemapper::prepareInterpKernelOnlyCC()
675 std::string srcMeth,trgMeth;
676 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
677 if(methodCpp!="P0P0")
678 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
679 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
680 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
681 const int srcMeshDim=src_mesh->getMeshDimension();
682 const int trgMeshDim=target_mesh->getMeshDimension();
683 if(trgMeshDim!=srcMeshDim)
684 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
689 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
690 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
691 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
692 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
697 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
698 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
699 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
700 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
705 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
706 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
707 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
708 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
712 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
714 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
716 _deno_multiply.clear();
717 _deno_multiply.resize(_matrix.size());
718 _deno_reverse_multiply.clear();
719 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
724 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss()
726 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
727 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 !");
728 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
729 const double *trgLocPtr=trgLoc->begin();
730 int trgSpaceDim=trgLoc->getNumberOfComponents();
731 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
732 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
734 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
735 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
736 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
737 throw INTERP_KERNEL::Exception(oss.str().c_str());
739 const int *srcOffsetArrPtr=srcOffsetArr->begin();
740 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
741 const double *srcLocPtr=srcLoc->begin();
742 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsArr,eltsIndexArr;
743 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
744 _matrix.resize(trgNbOfGaussPts);
745 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),eltsArr,eltsIndexArr);
746 const int *elts(eltsArr->begin()),*eltsIndex(eltsIndexArr->begin());
747 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfSrcCellsShTrgPts(eltsIndexArr->deltaShiftIndex());
748 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->getIdsNotEqual(0);
749 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
751 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
752 int srcCellId=elts[eltsIndex[*trgId]];
753 double dist=std::numeric_limits<double>::max();
755 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
757 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
759 for(int i=0;i<trgSpaceDim;i++)
760 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
762 { dist=tmp; srcEntry=srcId; }
764 _matrix[*trgId][srcEntry]=1.;
766 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
768 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->getIdsEqual(0);
769 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
770 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
771 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
772 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
773 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
775 _deno_multiply.clear();
776 _deno_multiply.resize(_matrix.size());
777 _deno_reverse_multiply.clear();
778 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
784 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
785 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
787 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method)
789 if(method=="GAUSSGAUSS")
791 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
792 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
793 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
794 throw INTERP_KERNEL::Exception(oss.str().c_str());
798 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
799 * 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
800 * only INTERP_KERNEL method should be applied.
802 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const
804 std::string srcm,trgm,method;
805 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
806 switch(_interp_matrix_pol)
808 case IK_ONLY_PREFERED:
812 std::string tmp1,tmp2;
813 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method.c_str(),tmp1,tmp2);
816 catch(INTERP_KERNEL::Exception& /*e*/)
821 case NOT_IK_ONLY_PREFERED:
825 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
828 catch(INTERP_KERNEL::Exception& /*e*/)
835 case NOT_IK_ONLY_FORCED:
838 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
842 void MEDCouplingRemapper::updateTime() const
846 void MEDCouplingRemapper::checkPrepare() const
848 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
850 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
851 if(!s->getMesh() || !t->getMesh())
852 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
856 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
857 * This method returns 3 informations (2 in ouput parameters and 1 in return).
859 * \param [out] srcMeth the string code of the discretization of source field template
860 * \param [out] trgMeth the string code of the discretization of target field template
861 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
863 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const
865 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
867 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
868 if(!s->getMesh() || !t->getMesh())
869 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
870 srcMeth=_src_ft->getDiscretization()->getRepr();
871 trgMeth=_target_ft->getDiscretization()->getRepr();
872 std::string method(srcMeth); method+=trgMeth;
876 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
880 if(matrixSuppression)
883 _deno_multiply.clear();
884 _deno_reverse_multiply.clear();
888 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue)
891 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
892 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
893 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
894 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
895 if(srcField->getNature()!=targetField->getNature())
896 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
897 DataArrayDouble *array=targetField->getArray();
898 int srcNbOfCompo=srcField->getNumberOfComponents();
901 if(srcNbOfCompo!=targetField->getNumberOfComponents())
902 throw INTERP_KERNEL::Exception("Number of components mismatch !");
907 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
908 array=DataArrayDouble::New();
909 array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
910 targetField->setArray(array);
913 computeDeno(srcField->getNature(),srcField,targetField);
914 double *resPointer=array->getPointer();
915 const double *inputPointer=srcField->getArray()->getConstPointer();
916 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
919 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
922 return computeDenoFromScratch(nat,srcField,trgField);
923 else if(nat!=_nature_of_deno)
924 return computeDenoFromScratch(nat,srcField,trgField);
925 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
926 return computeDenoFromScratch(nat,srcField,trgField);
929 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
932 _time_deno_update=getTimeOfThis();
933 switch(_nature_of_deno)
935 case ConservativeVolumic:
937 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
942 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
943 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
944 const double *denoPtr=deno->getArray()->getConstPointer();
945 const double *denoRPtr=denoR->getArray()->getConstPointer();
946 if(trgField->getMesh()->getMeshDimension()==-1)
948 double *denoRPtr2=denoR->getArray()->getPointer();
949 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
951 if(srcField->getMesh()->getMeshDimension()==-1)
953 double *denoPtr2=deno->getArray()->getPointer();
954 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
957 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
958 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
960 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
961 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
967 case IntegralGlobConstraint:
969 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
974 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
975 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
976 const double *denoPtr=deno->getArray()->getConstPointer();
977 const double *denoRPtr=denoR->getArray()->getConstPointer();
978 if(trgField->getMesh()->getMeshDimension()==-1)
980 double *denoRPtr2=denoR->getArray()->getPointer();
981 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
983 if(srcField->getMesh()->getMeshDimension()==-1)
985 double *denoPtr2=deno->getArray()->getPointer();
986 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
989 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
990 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
992 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
993 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
1000 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
1004 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
1007 double *tmp=new double[inputNbOfCompo];
1008 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1010 if((*iter1).empty())
1013 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1017 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1018 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1019 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1021 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1022 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1028 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1030 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1032 double *tmp=new double[inputNbOfCompo];
1033 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1034 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1036 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1038 isReached[(*iter2).first]=true;
1039 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1040 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1045 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1047 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1050 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1052 matOut.resize(nbColsMatIn);
1054 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1055 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1056 matOut[(*iter2).first][id]=(*iter2).second;
1059 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1060 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1062 std::map<int,double> values;
1064 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1067 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1069 sum+=(*iter2).second;
1070 values[(*iter2).first]+=(*iter2).second;
1072 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1073 deno[idx][(*iter2).first]=sum;
1076 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1078 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1079 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1083 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1084 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1086 std::map<int,double> values;
1088 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1091 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1093 sum+=(*iter2).second;
1094 values[(*iter2).first]+=(*iter2).second;
1096 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1097 denoReverse[(*iter2).first][idx]=sum;
1100 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1102 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1103 deno[idx][(*iter2).first]=values[(*iter2).first];
1107 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1108 const std::vector< std::map<int,double> >& m2D,
1109 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1110 const int *corrCellIdTrg)
1112 int nbOf2DCellsTrg=m2D.size();
1113 int nbOf1DCellsTrg=m1D.size();
1114 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1115 _matrix.resize(nbOf3DCellsTrg);
1117 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1119 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1122 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1124 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1126 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1133 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1136 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1138 std::cout << "Target Cell # " << id << " : ";
1139 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1140 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1141 std::cout << std::endl;
1145 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1151 * Returns the number of columns of matrix returned by MEDCouplingRemapper::getCrudeMatrix method.
1153 int MEDCouplingRemapper::getNumberOfColsOfMatrix() const
1155 return (int)_deno_reverse_multiply.size();
1159 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1160 * If not the behaviour is unpredictable.
1161 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1162 * 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.
1163 * 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
1166 * \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.
1167 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1168 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1170 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs)
1173 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1175 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1177 std::map<int,double>& rowNew=matrixNew[i];
1178 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1180 if(fabs((*it2).second)>maxValAbs)
1181 rowNew[(*it2).first]=(*it2).second;
1192 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1193 * If not the behaviour is unpredictable.
1194 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1195 * 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.
1196 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1197 * 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
1200 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1201 * \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
1202 * that all coefficients are null.
1203 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1205 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor)
1207 double maxVal=getMaxValueInCrudeMatrix();
1210 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1214 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1215 * If not the behaviour is unpredictable.
1216 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1217 * The returned value is positive.
1219 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const
1222 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1223 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1224 if(fabs((*it2).second)>ret)
1225 ret=fabs((*it2).second);