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) throw(INTERP_KERNEL::Exception)
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) throw(INTERP_KERNEL::Exception)
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() throw(INTERP_KERNEL::Exception)
82 int meshInterpType=((int)_src_ft->getMesh()->getType()*16)+(int)_target_ft->getMesh()->getType();
83 switch(meshInterpType)
85 case 85://Unstructured-Unstructured
86 return prepareInterpKernelOnlyUU();
87 case 87://Unstructured-Cartesian
88 return prepareInterpKernelOnlyUC();
89 case 117://Cartesian-Unstructured
90 return prepareInterpKernelOnlyCU();
91 case 119://Cartesian-Cartesian
92 return prepareInterpKernelOnlyCC();
93 case 136://Extruded-Extruded
94 return prepareInterpKernelOnlyEE();
96 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnly : Not managed type of meshes ! Dealt meshes type are : Unstructured<->Unstructured, Unstructured<->Cartesian, Cartesian<->Cartesian, Extruded<->Extruded !");
100 int MEDCouplingRemapper::prepareNotInterpKernelOnly() throw(INTERP_KERNEL::Exception)
102 std::string srcm,trgm,method;
103 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
104 switch(CheckInterpolationMethodManageableByNotOnlyInterpKernel(method))
107 return prepareNotInterpKernelOnlyGaussGauss();
110 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnly : INTERNAL ERROR ! the method \"" << method << "\" declared as managed bu not implemented !";
111 throw INTERP_KERNEL::Exception(oss.str().c_str());
117 * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
118 * 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.
120 * \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.
121 * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
123 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
125 transferUnderground(srcField,targetField,true,dftValue);
129 * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
130 * 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
132 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
134 * \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.
135 * \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.
137 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
139 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
142 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
145 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
146 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
147 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
148 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
149 if(srcField->getNature()!=targetField->getNature())
150 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
151 DataArrayDouble *array=srcField->getArray();
152 int trgNbOfCompo=targetField->getNumberOfComponents();
155 if(trgNbOfCompo!=srcField->getNumberOfComponents())
156 throw INTERP_KERNEL::Exception("Number of components mismatch !");
160 array=DataArrayDouble::New();
161 array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
162 srcField->setArray(array);
165 computeDeno(srcField->getNature(),srcField,targetField);
166 double *resPointer=array->getPointer();
167 const double *inputPointer=targetField->getArray()->getConstPointer();
168 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
171 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
174 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
175 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
176 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_target_ft,srcField->getTimeDiscretization());
177 ret->setNature(srcField->getNature());
178 transfer(srcField,ret,dftValue);
179 ret->copyAllTinyAttrFrom(srcField);//perform copy of tiny strings after and not before transfer because the array will be created on transfer
183 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
186 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
187 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
188 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(*_src_ft,targetField->getTimeDiscretization());
189 ret->setNature(targetField->getNature());
190 reverseTransfer(ret,targetField,dftValue);
191 ret->copyAllTinyAttrFrom(targetField);//perform copy of tiny strings after and not before reverseTransfer because the array will be created on reverseTransfer
196 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
197 * is here only for automatic CORBA generators.
199 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
201 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
205 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
206 * is here only for automatic CORBA generators.
208 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
210 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
214 * This method does nothing more than inherited INTERP_KERNEL::InterpolationOptions::setOptionInt method. This method
215 * is here only for automatic CORBA generators.
217 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
219 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
223 * This method returns the interpolation matrix policy. This policy specifies which interpolation matrix method to keep or prefered.
224 * If interpolation matrix policy is :
226 * - 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
227 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
228 * If not, the \b not only INTERP_KERNEL method will be attempt.
230 * - 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
231 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
232 * If not, the INTERP_KERNEL only method will be attempt.
234 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
236 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
238 * \sa MEDCouplingRemapper::setInterpolationMatrixPolicy
240 int MEDCouplingRemapper::getInterpolationMatrixPolicy() const
242 return _interp_matrix_pol;
246 * 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
247 * 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
248 * [0 (IK_PREFERED) , 1 (NOT_IK_PREFERED), 2 (IK_ONLY_FORCED), 3 (NOT_IK_ONLY_FORCED)].
250 * If interpolation matrix policy is :
252 * - 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
253 * regarding spatial discretization of source and target with INTERP_KERNEL only method, INTERP_KERNEL only method will be performed.
254 * If not, the \b not only INTERP_KERNEL method will be attempt.
256 * - 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
257 * regarding spatial discretization of source and target with \b NOT only INTERP_KERNEL method, \b NOT only INTERP_KERNEL method, will be performed.
258 * If not, the INTERP_KERNEL only method will be attempt.
260 * - IK_ONLY_FORCED (2) : Only INTERP_KERNEL only method will be launched.
262 * - NOT_IK_ONLY_FORCED (3) : Only \b NOT INTERP_KERNEL only method will be launched.
264 * \input newInterpMatPol the new interpolation matrix method policy. This parameter is of type \c int and not of type \c ParaMEDMEM::InterpolationMatrixPolicy
265 * for automatic generation of CORBA component.
267 * \sa MEDCouplingRemapper::getInterpolationMatrixPolicy
269 void MEDCouplingRemapper::setInterpolationMatrixPolicy(int newInterpMatPol) throw(INTERP_KERNEL::Exception)
271 switch(newInterpMatPol)
274 _interp_matrix_pol=IK_ONLY_PREFERED;
277 _interp_matrix_pol=NOT_IK_ONLY_PREFERED;
280 _interp_matrix_pol=IK_ONLY_FORCED;
283 _interp_matrix_pol=NOT_IK_ONLY_FORCED;
286 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 !");
290 int MEDCouplingRemapper::prepareInterpKernelOnlyUU() throw(INTERP_KERNEL::Exception)
292 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
293 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
294 std::string srcMeth,trgMeth;
295 std::string method=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
296 const int srcMeshDim=src_mesh->getMeshDimension();
299 srcSpaceDim=src_mesh->getSpaceDimension();
300 const int trgMeshDim=target_mesh->getMeshDimension();
303 trgSpaceDim=target_mesh->getSpaceDimension();
304 if(trgSpaceDim!=srcSpaceDim)
305 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
306 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
308 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
310 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
311 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
312 INTERP_KERNEL::Interpolation1D interpolation(*this);
313 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
315 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
317 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
318 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
319 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
320 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
322 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
324 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
325 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
326 INTERP_KERNEL::Interpolation2D interpolation(*this);
327 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
329 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
331 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
332 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
333 INTERP_KERNEL::Interpolation3D interpolation(*this);
334 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
336 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
338 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
339 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
340 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
341 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
343 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
345 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
346 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
347 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
348 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
349 INTERP_KERNEL::Interpolation3D interpolation(*this);
350 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
352 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
354 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
355 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
356 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
357 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
358 INTERP_KERNEL::Interpolation3D interpolation(*this);
359 std::vector<std::map<int,double> > matrixTmp;
360 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
361 ReverseMatrix(matrixTmp,nbCols,_matrix);
362 nbCols=matrixTmp.size();
364 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
366 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
368 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
369 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
370 INTERP_KERNEL::Interpolation2D interpolation(*this);
371 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
375 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
376 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
377 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
378 std::vector<std::map<int,double> > matrixTmp;
379 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
380 ReverseMatrix(matrixTmp,nbCols,_matrix);
381 nbCols=matrixTmp.size();
382 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
383 if(!duplicateFaces.empty())
385 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
386 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
388 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
389 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
395 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
397 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
399 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
400 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
401 INTERP_KERNEL::Interpolation2D interpolation(*this);
402 std::vector<std::map<int,double> > matrixTmp;
403 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
404 ReverseMatrix(matrixTmp,nbCols,_matrix);
405 nbCols=matrixTmp.size();
409 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
410 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
411 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
412 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
413 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
414 if(!duplicateFaces.empty())
416 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
417 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
419 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
420 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
426 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
428 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
429 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
430 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
431 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method.c_str());
432 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
433 if(!duplicateFaces.empty())
435 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
436 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
438 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
439 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
444 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
446 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
447 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
448 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
449 std::vector<std::map<int,double> > matrixTmp;
450 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method.c_str());
451 ReverseMatrix(matrixTmp,nbCols,_matrix);
452 nbCols=matrixTmp.size();
453 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
454 if(!duplicateFaces.empty())
456 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
457 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
459 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
460 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
465 else if(trgMeshDim==-1)
467 if(srcMeshDim==2 && srcSpaceDim==2)
469 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
470 INTERP_KERNEL::Interpolation2D interpolation(*this);
471 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
473 else if(srcMeshDim==3 && srcSpaceDim==3)
475 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
476 INTERP_KERNEL::Interpolation3D interpolation(*this);
477 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
479 else if(srcMeshDim==2 && srcSpaceDim==3)
481 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
482 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
483 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,srcMeth.c_str());
486 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
488 else if(srcMeshDim==-1)
490 if(trgMeshDim==2 && trgSpaceDim==2)
492 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
493 INTERP_KERNEL::Interpolation2D interpolation(*this);
494 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
496 else if(trgMeshDim==3 && trgSpaceDim==3)
498 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
499 INTERP_KERNEL::Interpolation3D interpolation(*this);
500 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
502 else if(trgMeshDim==2 && trgSpaceDim==3)
504 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
505 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
506 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,trgMeth.c_str());
509 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
512 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
513 _deno_multiply.clear();
514 _deno_multiply.resize(_matrix.size());
515 _deno_reverse_multiply.clear();
516 _deno_reverse_multiply.resize(nbCols);
521 int MEDCouplingRemapper::prepareInterpKernelOnlyEE() throw(INTERP_KERNEL::Exception)
523 std::string srcMeth,trgMeth;
524 std::string methC=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
525 const MEDCouplingExtrudedMesh *src_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_src_ft->getMesh());
526 const MEDCouplingExtrudedMesh *target_mesh=static_cast<const MEDCouplingExtrudedMesh *>(_target_ft->getMesh());
528 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyEE : Only P0P0 method implemented for Extruded/Extruded meshes !");
529 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
530 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
531 INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
532 std::vector<std::map<int,double> > matrix2D;
533 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,methC.c_str());
534 MEDCouplingUMesh *s1D,*t1D;
536 MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
537 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
538 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
539 std::vector<std::map<int,double> > matrix1D;
540 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
541 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,methC.c_str());
544 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
545 target_mesh->getMesh3DIds()->getConstPointer());
547 _deno_multiply.clear();
548 _deno_multiply.resize(_matrix.size());
549 _deno_reverse_multiply.clear();
550 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
555 int MEDCouplingRemapper::prepareInterpKernelOnlyUC() throw(INTERP_KERNEL::Exception)
557 std::string srcMeth,trgMeth;
558 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
559 if(methodCpp!="P0P0")
560 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only P0P0 interpolation supported for the moment !");
561 const MEDCouplingUMesh *src_mesh=static_cast<const MEDCouplingUMesh *>(_src_ft->getMesh());
562 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
563 const int srcMeshDim=src_mesh->getMeshDimension();
564 const int srcSpceDim=src_mesh->getSpaceDimension();
565 const int trgMeshDim=target_mesh->getMeshDimension();
566 if(srcMeshDim!=srcSpceDim || srcMeshDim!=trgMeshDim)
567 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 !");
568 std::vector<std::map<int,double> > res;
573 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
574 MEDCouplingNormalizedUnstructuredMesh<1,1> sourceWrapper(src_mesh);
575 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
576 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
581 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
582 MEDCouplingNormalizedUnstructuredMesh<2,2> sourceWrapper(src_mesh);
583 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
584 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
589 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
590 MEDCouplingNormalizedUnstructuredMesh<3,3> sourceWrapper(src_mesh);
591 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
592 myInterpolator.interpolateMeshes(targetWrapper,sourceWrapper,res,"P0P0");
596 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyUC : only dimension 1 2 or 3 supported !");
598 ReverseMatrix(res,target_mesh->getNumberOfCells(),_matrix);
599 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
601 _deno_multiply.clear();
602 _deno_multiply.resize(_matrix.size());
603 _deno_reverse_multiply.clear();
604 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
609 int MEDCouplingRemapper::prepareInterpKernelOnlyCU() throw(INTERP_KERNEL::Exception)
611 std::string srcMeth,trgMeth;
612 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
613 if(methodCpp!="P0P0")
614 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only P0P0 interpolation supported for the moment !");
615 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
616 const MEDCouplingUMesh *target_mesh=static_cast<const MEDCouplingUMesh *>(_target_ft->getMesh());
617 const int srcMeshDim=src_mesh->getMeshDimension();
618 const int trgMeshDim=target_mesh->getMeshDimension();
619 const int trgSpceDim=target_mesh->getSpaceDimension();
620 if(trgMeshDim!=trgSpceDim || trgMeshDim!=srcMeshDim)
621 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 !");
626 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
627 MEDCouplingNormalizedUnstructuredMesh<1,1> targetWrapper(target_mesh);
628 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
629 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
634 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
635 MEDCouplingNormalizedUnstructuredMesh<2,2> targetWrapper(target_mesh);
636 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
637 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
642 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
643 MEDCouplingNormalizedUnstructuredMesh<3,3> targetWrapper(target_mesh);
644 INTERP_KERNEL::InterpolationCU myInterpolator(*this);
645 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
649 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCU : only dimension 1 2 or 3 supported !");
651 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
653 _deno_multiply.clear();
654 _deno_multiply.resize(_matrix.size());
655 _deno_reverse_multiply.clear();
656 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
661 int MEDCouplingRemapper::prepareInterpKernelOnlyCC() throw(INTERP_KERNEL::Exception)
663 std::string srcMeth,trgMeth;
664 std::string methodCpp=checkAndGiveInterpolationMethodStr(srcMeth,trgMeth);
665 if(methodCpp!="P0P0")
666 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only P0P0 interpolation supported for the moment !");
667 const MEDCouplingCMesh *src_mesh=static_cast<const MEDCouplingCMesh *>(_src_ft->getMesh());
668 const MEDCouplingCMesh *target_mesh=static_cast<const MEDCouplingCMesh *>(_target_ft->getMesh());
669 const int srcMeshDim=src_mesh->getMeshDimension();
670 const int trgMeshDim=target_mesh->getMeshDimension();
671 if(trgMeshDim!=srcMeshDim)
672 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : dim of target cartesian should be equal to dim of source cartesian dimension !");
677 MEDCouplingNormalizedCartesianMesh<1> sourceWrapper(src_mesh);
678 MEDCouplingNormalizedCartesianMesh<1> targetWrapper(target_mesh);
679 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
680 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
685 MEDCouplingNormalizedCartesianMesh<2> sourceWrapper(src_mesh);
686 MEDCouplingNormalizedCartesianMesh<2> targetWrapper(target_mesh);
687 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
688 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
693 MEDCouplingNormalizedCartesianMesh<3> sourceWrapper(src_mesh);
694 MEDCouplingNormalizedCartesianMesh<3> targetWrapper(target_mesh);
695 INTERP_KERNEL::InterpolationCC myInterpolator(*this);
696 myInterpolator.interpolateMeshes(sourceWrapper,targetWrapper,_matrix,"P0P0");
700 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::prepareInterpKernelOnlyCC : only dimension 1 2 or 3 supported !");
702 nullifiedTinyCoeffInCrudeMatrixAbs(0.);
704 _deno_multiply.clear();
705 _deno_multiply.resize(_matrix.size());
706 _deno_reverse_multiply.clear();
707 _deno_reverse_multiply.resize(src_mesh->getNumberOfCells());
712 int MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss() throw(INTERP_KERNEL::Exception)
714 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
715 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 !");
716 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> trgLoc=_target_ft->getLocalizationOfDiscr();
717 const double *trgLocPtr=trgLoc->begin();
718 int trgSpaceDim=trgLoc->getNumberOfComponents();
719 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcOffsetArr=_src_ft->getDiscretization()->getOffsetArr(_src_ft->getMesh());
720 if(trgSpaceDim!=_src_ft->getMesh()->getSpaceDimension())
722 std::ostringstream oss; oss << "MEDCouplingRemapper::prepareNotInterpKernelOnlyGaussGauss : space dimensions mismatch between source and target !";
723 oss << " Target discretization localization has dimension " << trgSpaceDim << ", whereas the space dimension of source is equal to ";
724 oss << _src_ft->getMesh()->getSpaceDimension() << " !";
725 throw INTERP_KERNEL::Exception(oss.str().c_str());
727 const int *srcOffsetArrPtr=srcOffsetArr->begin();
728 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> srcLoc=_src_ft->getLocalizationOfDiscr();
729 const double *srcLocPtr=srcLoc->begin();
730 std::vector<int> elts,eltsIndex;
731 int trgNbOfGaussPts=trgLoc->getNumberOfTuples();
732 _matrix.resize(trgNbOfGaussPts);
733 _src_ft->getMesh()->getCellsContainingPoints(trgLoc->begin(),trgNbOfGaussPts,getPrecision(),elts,eltsIndex);
734 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> eltsIndex2=DataArrayInt::New(); eltsIndex2->useArray(&eltsIndex[0],false,CPP_DEALLOC,(int)eltsIndex.size(),1);
735 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> nbOfSrcCellsShTrgPts=eltsIndex2->deltaShiftIndex();
736 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> ids0=nbOfSrcCellsShTrgPts->getIdsNotEqual(0);
737 for(const int *trgId=ids0->begin();trgId!=ids0->end();trgId++)
739 const double *ptTrg=trgLocPtr+trgSpaceDim*(*trgId);
740 int srcCellId=elts[eltsIndex[*trgId]];
741 double dist=std::numeric_limits<double>::max();
743 for(int srcId=srcOffsetArrPtr[srcCellId];srcId<srcOffsetArrPtr[srcCellId+1];srcId++)
745 const double *ptSrc=srcLocPtr+trgSpaceDim*srcId;
747 for(int i=0;i<trgSpaceDim;i++)
748 tmp+=(ptTrg[i]-ptSrc[i])*(ptTrg[i]-ptSrc[i]);
750 { dist=tmp; srcEntry=srcId; }
752 _matrix[*trgId][srcEntry]=1.;
754 if(ids0->getNumberOfTuples()!=trgNbOfGaussPts)
756 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> orphanTrgIds=nbOfSrcCellsShTrgPts->getIdsEqual(0);
757 MEDCouplingAutoRefCountObjectPtr<DataArrayDouble> orphanTrg=trgLoc->selectByTupleId(orphanTrgIds->begin(),orphanTrgIds->end());
758 MEDCouplingAutoRefCountObjectPtr<DataArrayInt> srcIdPerTrg=srcLoc->findClosestTupleId(orphanTrg);
759 const int *srcIdPerTrgPtr=srcIdPerTrg->begin();
760 for(const int *orphanTrgId=orphanTrgIds->begin();orphanTrgId!=orphanTrgIds->end();orphanTrgId++,srcIdPerTrgPtr++)
761 _matrix[*orphanTrgId][*srcIdPerTrgPtr]=2.;
763 _deno_multiply.clear();
764 _deno_multiply.resize(_matrix.size());
765 _deno_reverse_multiply.clear();
766 _deno_reverse_multiply.resize(srcLoc->getNumberOfTuples());
772 * This method checks that the input interpolation \a method is managed by not INTERP_KERNEL only methods.
773 * If no an INTERP_KERNEL::Exception will be thrown. If yes, a magic number will be returned to switch in the MEDCouplingRemapper::prepareNotInterpKernelOnly method.
775 int MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel(const std::string& method) throw(INTERP_KERNEL::Exception)
777 if(method=="GAUSSGAUSS")
779 std::ostringstream oss; oss << "MEDCouplingRemapper::CheckInterpolationMethodManageableByNotOnlyInterpKernel : ";
780 oss << "The method \"" << method << "\" is not manageable by not INTERP_KERNEL only method.";
781 oss << " Not only INTERP_KERNEL methods dealed are : GAUSSGAUSS !";
782 throw INTERP_KERNEL::Exception(oss.str().c_str());
786 * This method determines regarding \c _interp_matrix_pol attribute ( set by MEDCouplingRemapper::setInterpolationMatrixPolicy and by default equal
787 * 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
788 * only INTERP_KERNEL method should be applied.
790 bool MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly() const throw(INTERP_KERNEL::Exception)
792 std::string srcm,trgm,method;
793 method=checkAndGiveInterpolationMethodStr(srcm,trgm);
794 switch(_interp_matrix_pol)
796 case IK_ONLY_PREFERED:
800 std::string tmp1,tmp2;
801 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::CheckAndSplitInterpolationMethod(method.c_str(),tmp1,tmp2);
804 catch(INTERP_KERNEL::Exception& e)
809 case NOT_IK_ONLY_PREFERED:
813 CheckInterpolationMethodManageableByNotOnlyInterpKernel(method);
816 catch(INTERP_KERNEL::Exception& e)
823 case NOT_IK_ONLY_FORCED:
826 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::isInterpKernelOnlyOrNotOnly : internal error ! The interpolation matrix policy is not managed ! Try to change it using MEDCouplingRemapper::setInterpolationMatrixPolicy !");
830 void MEDCouplingRemapper::updateTime() const
834 void MEDCouplingRemapper::checkPrepare() const throw(INTERP_KERNEL::Exception)
836 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
838 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that MEDCouplingRemapper::prepare(Ex) has not been called !");
839 if(!s->getMesh() || !t->getMesh())
840 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkPrepare : it appears that no all field templates have their mesh set !");
844 * This method builds a code considering already set field discretization int \a this : \a _src_ft and \a _target_ft.
845 * This method returns 3 informations (2 in ouput parameters and 1 in return).
847 * \param [out] srcMeth the string code of the discretization of source field template
848 * \param [out] trgMeth the string code of the discretization of target field template
849 * \return the standardized string code (compatible with INTERP_KERNEL) for matrix of numerators (in \a _matrix)
851 std::string MEDCouplingRemapper::checkAndGiveInterpolationMethodStr(std::string& srcMeth, std::string& trgMeth) const throw(INTERP_KERNEL::Exception)
853 const MEDCouplingFieldTemplate *s(_src_ft),*t(_target_ft);
855 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have been set !");
856 if(!s->getMesh() || !t->getMesh())
857 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::checkAndGiveInterpolationMethodStr : it appears that no all field templates have their mesh set !");
858 srcMeth=_src_ft->getDiscretization()->getRepr();
859 trgMeth=_target_ft->getDiscretization()->getRepr();
860 std::string method(srcMeth); method+=trgMeth;
864 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
868 if(matrixSuppression)
871 _deno_multiply.clear();
872 _deno_reverse_multiply.clear();
876 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
879 if(_src_ft->getDiscretization()->getStringRepr()!=srcField->getDiscretization()->getStringRepr())
880 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
881 if(_target_ft->getDiscretization()->getStringRepr()!=targetField->getDiscretization()->getStringRepr())
882 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
883 if(srcField->getNature()!=targetField->getNature())
884 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
885 DataArrayDouble *array=targetField->getArray();
886 int srcNbOfCompo=srcField->getNumberOfComponents();
889 if(srcNbOfCompo!=targetField->getNumberOfComponents())
890 throw INTERP_KERNEL::Exception("Number of components mismatch !");
895 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
896 array=DataArrayDouble::New();
897 array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
898 targetField->setArray(array);
901 computeDeno(srcField->getNature(),srcField,targetField);
902 double *resPointer=array->getPointer();
903 const double *inputPointer=srcField->getArray()->getConstPointer();
904 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
907 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
910 return computeDenoFromScratch(nat,srcField,trgField);
911 else if(nat!=_nature_of_deno)
912 return computeDenoFromScratch(nat,srcField,trgField);
913 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
914 return computeDenoFromScratch(nat,srcField,trgField);
917 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
920 _time_deno_update=getTimeOfThis();
921 switch(_nature_of_deno)
923 case ConservativeVolumic:
925 ComputeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
930 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
931 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
932 const double *denoPtr=deno->getArray()->getConstPointer();
933 const double *denoRPtr=denoR->getArray()->getConstPointer();
934 if(trgField->getMesh()->getMeshDimension()==-1)
936 double *denoRPtr2=denoR->getArray()->getPointer();
937 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
939 if(srcField->getMesh()->getMeshDimension()==-1)
941 double *denoPtr2=deno->getArray()->getPointer();
942 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
945 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
946 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
948 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
949 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
955 case IntegralGlobConstraint:
957 ComputeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
962 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),getMeasureAbsStatus());
963 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),getMeasureAbsStatus());
964 const double *denoPtr=deno->getArray()->getConstPointer();
965 const double *denoRPtr=denoR->getArray()->getConstPointer();
966 if(trgField->getMesh()->getMeshDimension()==-1)
968 double *denoRPtr2=denoR->getArray()->getPointer();
969 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
971 if(srcField->getMesh()->getMeshDimension()==-1)
973 double *denoPtr2=deno->getArray()->getPointer();
974 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
977 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
978 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
980 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
981 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
988 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
992 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
995 double *tmp=new double[inputNbOfCompo];
996 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1001 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1005 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
1006 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
1007 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
1009 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
1010 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
1016 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
1018 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
1020 double *tmp=new double[inputNbOfCompo];
1021 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
1022 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
1024 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1026 isReached[(*iter2).first]=true;
1027 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
1028 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
1033 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
1035 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
1038 void MEDCouplingRemapper::ReverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
1040 matOut.resize(nbColsMatIn);
1042 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
1043 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1044 matOut[(*iter2).first][id]=(*iter2).second;
1047 void MEDCouplingRemapper::ComputeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
1048 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1050 std::map<int,double> values;
1052 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1055 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1057 sum+=(*iter2).second;
1058 values[(*iter2).first]+=(*iter2).second;
1060 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1061 deno[idx][(*iter2).first]=sum;
1064 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1066 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1067 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
1071 void MEDCouplingRemapper::ComputeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
1072 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
1074 std::map<int,double> values;
1076 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1079 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1081 sum+=(*iter2).second;
1082 values[(*iter2).first]+=(*iter2).second;
1084 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1085 denoReverse[(*iter2).first][idx]=sum;
1088 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
1090 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1091 deno[idx][(*iter2).first]=values[(*iter2).first];
1095 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
1096 const std::vector< std::map<int,double> >& m2D,
1097 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
1098 const int *corrCellIdTrg)
1100 int nbOf2DCellsTrg=m2D.size();
1101 int nbOf1DCellsTrg=m1D.size();
1102 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
1103 _matrix.resize(nbOf3DCellsTrg);
1105 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
1107 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
1110 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
1112 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
1114 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
1121 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
1124 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
1126 std::cout << "Target Cell # " << id << " : ";
1127 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
1128 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
1129 std::cout << std::endl;
1133 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
1139 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1140 * If not the behaviour is unpredictable.
1141 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
1142 * 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.
1143 * 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
1146 * \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.
1147 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
1148 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
1150 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception)
1153 std::vector<std::map<int,double> > matrixNew(_matrix.size());
1155 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
1157 std::map<int,double>& rowNew=matrixNew[i];
1158 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1160 if(fabs((*it2).second)>maxValAbs)
1161 rowNew[(*it2).first]=(*it2).second;
1172 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1173 * If not the behaviour is unpredictable.
1174 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
1175 * 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.
1176 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
1177 * 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
1180 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
1181 * \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
1182 * that all coefficients are null.
1183 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
1185 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception)
1187 double maxVal=getMaxValueInCrudeMatrix();
1190 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
1194 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
1195 * If not the behaviour is unpredictable.
1196 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
1197 * The returned value is positive.
1199 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception)
1202 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
1203 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
1204 if(fabs((*it2).second)>ret)
1205 ret=fabs((*it2).second);