1 // Copyright (C) 2007-2012 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 "MEDCouplingNormalizedUnstructuredMesh.txx"
29 #include "Interpolation1D.txx"
30 #include "Interpolation2DCurve.hxx"
31 #include "Interpolation2D.txx"
32 #include "Interpolation3D.txx"
33 #include "Interpolation3DSurf.hxx"
34 #include "Interpolation2D1D.txx"
35 #include "Interpolation3D2D.txx"
37 using namespace ParaMEDMEM;
39 MEDCouplingRemapper::MEDCouplingRemapper():_src_mesh(0),_target_mesh(0),_nature_of_deno(NoNature),_time_deno_update(0)
43 MEDCouplingRemapper::~MEDCouplingRemapper()
48 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const char *method) throw(INTERP_KERNEL::Exception)
51 _src_mesh=const_cast<MEDCouplingMesh *>(srcMesh); _target_mesh=const_cast<MEDCouplingMesh *>(targetMesh);
52 _src_mesh->incrRef(); _target_mesh->incrRef();
53 int meshInterpType=((int)_src_mesh->getType()*16)+(int)_target_mesh->getType();
54 switch(meshInterpType)
56 case 85://Unstructured-Unstructured
57 return prepareUU(method);
58 case 136://Extruded-Extruded
59 return prepareEE(method);
61 throw INTERP_KERNEL::Exception("Not managed type of meshes !");
65 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target) throw(INTERP_KERNEL::Exception)
67 std::string meth(src->getDiscretization()->getStringRepr());
68 meth+=target->getDiscretization()->getStringRepr();
69 return prepare(src->getMesh(),target->getMesh(),meth.c_str());
73 * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
74 * 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.
76 * \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.
77 * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
79 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
81 transferUnderground(srcField,targetField,true,dftValue);
85 * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
86 * 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
88 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
90 * \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.
91 * \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.
93 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
95 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
98 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
100 if(_src_method!=srcField->getDiscretization()->getStringRepr())
101 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
102 if(_target_method!=targetField->getDiscretization()->getStringRepr())
103 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
104 if(srcField->getNature()!=targetField->getNature())
105 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
106 DataArrayDouble *array=srcField->getArray();
107 int trgNbOfCompo=targetField->getNumberOfComponents();
110 if(trgNbOfCompo!=srcField->getNumberOfComponents())
111 throw INTERP_KERNEL::Exception("Number of components mismatch !");
115 array=DataArrayDouble::New();
116 array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
117 srcField->setArray(array);
120 computeDeno(srcField->getNature(),srcField,targetField);
121 double *resPointer=array->getPointer();
122 const double *inputPointer=targetField->getArray()->getConstPointer();
123 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
126 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
128 if(_src_method!=srcField->getDiscretization()->getStringRepr())
129 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
130 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),srcField->getTimeDiscretization());
131 ret->copyTinyAttrFrom(srcField);
132 ret->setNature(srcField->getNature());
133 ret->setMesh(_target_mesh);
134 transfer(srcField,ret,dftValue);
138 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
140 if(_target_method!=targetField->getDiscretization()->getStringRepr())
141 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
142 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),targetField->getTimeDiscretization());
143 ret->copyTinyAttrFrom(targetField);
144 ret->setNature(targetField->getNature());
145 ret->setMesh(_src_mesh);
146 reverseTransfer(ret,targetField,dftValue);
150 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
152 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
155 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
157 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
160 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
162 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
165 int MEDCouplingRemapper::prepareUU(const char *method) throw(INTERP_KERNEL::Exception)
167 MEDCouplingUMesh *src_mesh=(MEDCouplingUMesh *)_src_mesh;
168 MEDCouplingUMesh *target_mesh=(MEDCouplingUMesh *)_target_mesh;
169 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
170 const int srcMeshDim=src_mesh->getMeshDimension();
173 srcSpaceDim=src_mesh->getSpaceDimension();
174 const int trgMeshDim=target_mesh->getMeshDimension();
177 trgSpaceDim=target_mesh->getSpaceDimension();
178 if(trgSpaceDim!=srcSpaceDim)
179 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
180 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
182 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
184 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
185 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
186 INTERP_KERNEL::Interpolation1D interpolation(*this);
187 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
189 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
191 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
192 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
193 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
194 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
196 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
198 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
199 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
200 INTERP_KERNEL::Interpolation2D interpolation(*this);
201 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
203 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
205 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
206 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
207 INTERP_KERNEL::Interpolation3D interpolation(*this);
208 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
210 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
212 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
213 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
214 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
215 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
217 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
219 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
220 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
221 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
222 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
223 INTERP_KERNEL::Interpolation3D interpolation(*this);
224 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
226 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
228 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
229 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
230 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
231 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
232 INTERP_KERNEL::Interpolation3D interpolation(*this);
233 std::vector<std::map<int,double> > matrixTmp;
234 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
235 reverseMatrix(matrixTmp,nbCols,_matrix);
236 nbCols=matrixTmp.size();
238 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
240 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
242 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
243 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
244 INTERP_KERNEL::Interpolation2D interpolation(*this);
245 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
249 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
250 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
251 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
252 std::vector<std::map<int,double> > matrixTmp;
253 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
254 reverseMatrix(matrixTmp,nbCols,_matrix);
255 nbCols=matrixTmp.size();
256 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
257 if(!duplicateFaces.empty())
259 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
260 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
262 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
263 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
269 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
271 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
273 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
274 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
275 INTERP_KERNEL::Interpolation2D interpolation(*this);
276 std::vector<std::map<int,double> > matrixTmp;
277 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
278 reverseMatrix(matrixTmp,nbCols,_matrix);
279 nbCols=matrixTmp.size();
283 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
284 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
285 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
286 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
287 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
288 if(!duplicateFaces.empty())
290 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
291 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
293 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
294 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
300 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
302 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
303 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
304 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
305 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
306 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
307 if(!duplicateFaces.empty())
309 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
310 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
312 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
313 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
318 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
320 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
321 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
322 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
323 std::vector<std::map<int,double> > matrixTmp;
324 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
325 reverseMatrix(matrixTmp,nbCols,_matrix);
326 nbCols=matrixTmp.size();
327 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
328 if(!duplicateFaces.empty())
330 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
331 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
333 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
334 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
339 else if(trgMeshDim==-1)
341 if(srcMeshDim==2 && srcSpaceDim==2)
343 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
344 INTERP_KERNEL::Interpolation2D interpolation(*this);
345 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
347 else if(srcMeshDim==3 && srcSpaceDim==3)
349 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
350 INTERP_KERNEL::Interpolation3D interpolation(*this);
351 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
353 else if(srcMeshDim==2 && srcSpaceDim==3)
355 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
356 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
357 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
360 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
362 else if(srcMeshDim==-1)
364 if(trgMeshDim==2 && trgSpaceDim==2)
366 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
367 INTERP_KERNEL::Interpolation2D interpolation(*this);
368 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
370 else if(trgMeshDim==3 && trgSpaceDim==3)
372 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
373 INTERP_KERNEL::Interpolation3D interpolation(*this);
374 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
376 else if(trgMeshDim==2 && trgSpaceDim==3)
378 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
379 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
380 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
383 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
386 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
387 _deno_multiply.clear();
388 _deno_multiply.resize(_matrix.size());
389 _deno_reverse_multiply.clear();
390 _deno_reverse_multiply.resize(nbCols);
395 int MEDCouplingRemapper::prepareEE(const char *method) throw(INTERP_KERNEL::Exception)
397 MEDCouplingExtrudedMesh *src_mesh=(MEDCouplingExtrudedMesh *)_src_mesh;
398 MEDCouplingExtrudedMesh *target_mesh=(MEDCouplingExtrudedMesh *)_target_mesh;
399 std::string methC(method);
401 throw INTERP_KERNEL::Exception("Only P0P0 method implemented for Extruded/Extruded meshes !");
402 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
403 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
404 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
405 INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
406 std::vector<std::map<int,double> > matrix2D;
407 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
408 MEDCouplingUMesh *s1D,*t1D;
410 MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
411 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
412 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
413 std::vector<std::map<int,double> > matrix1D;
414 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
415 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
418 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
419 target_mesh->getMesh3DIds()->getConstPointer());
421 _deno_multiply.clear();
422 _deno_multiply.resize(_matrix.size());
423 _deno_reverse_multiply.clear();
424 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
429 void MEDCouplingRemapper::updateTime() const
433 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
436 _src_mesh->decrRef();
438 _target_mesh->decrRef();
441 if(matrixSuppression)
444 _deno_multiply.clear();
445 _deno_reverse_multiply.clear();
449 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
451 if(_src_method!=srcField->getDiscretization()->getStringRepr())
452 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
453 if(_target_method!=targetField->getDiscretization()->getStringRepr())
454 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
455 if(srcField->getNature()!=targetField->getNature())
456 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
457 DataArrayDouble *array=targetField->getArray();
458 int srcNbOfCompo=srcField->getNumberOfComponents();
461 if(srcNbOfCompo!=targetField->getNumberOfComponents())
462 throw INTERP_KERNEL::Exception("Number of components mismatch !");
467 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
468 array=DataArrayDouble::New();
469 array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
470 targetField->setArray(array);
473 computeDeno(srcField->getNature(),srcField,targetField);
474 double *resPointer=array->getPointer();
475 const double *inputPointer=srcField->getArray()->getConstPointer();
476 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
479 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
482 return computeDenoFromScratch(nat,srcField,trgField);
483 else if(nat!=_nature_of_deno)
484 return computeDenoFromScratch(nat,srcField,trgField);
485 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
486 return computeDenoFromScratch(nat,srcField,trgField);
489 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
492 _time_deno_update=getTimeOfThis();
493 switch(_nature_of_deno)
495 case ConservativeVolumic:
497 computeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
502 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
503 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
504 const double *denoPtr=deno->getArray()->getConstPointer();
505 const double *denoRPtr=denoR->getArray()->getConstPointer();
506 if(trgField->getMesh()->getMeshDimension()==-1)
508 double *denoRPtr2=denoR->getArray()->getPointer();
509 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
511 if(srcField->getMesh()->getMeshDimension()==-1)
513 double *denoPtr2=deno->getArray()->getPointer();
514 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
517 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
518 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
520 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
521 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
527 case IntegralGlobConstraint:
529 computeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
534 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
535 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
536 const double *denoPtr=deno->getArray()->getConstPointer();
537 const double *denoRPtr=denoR->getArray()->getConstPointer();
538 if(trgField->getMesh()->getMeshDimension()==-1)
540 double *denoRPtr2=denoR->getArray()->getPointer();
541 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
543 if(srcField->getMesh()->getMeshDimension()==-1)
545 double *denoPtr2=deno->getArray()->getPointer();
546 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
549 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
550 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
552 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
553 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
560 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
564 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
567 double *tmp=new double[inputNbOfCompo];
568 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
573 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
577 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
578 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
579 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
581 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
582 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
588 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
590 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
592 double *tmp=new double[inputNbOfCompo];
593 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
594 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
596 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
598 isReached[(*iter2).first]=true;
599 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
600 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
605 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
607 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
610 void MEDCouplingRemapper::reverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
612 matOut.resize(nbColsMatIn);
614 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
615 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
616 matOut[(*iter2).first][id]=(*iter2).second;
619 void MEDCouplingRemapper::computeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
620 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
622 std::map<int,double> values;
624 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
627 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
629 sum+=(*iter2).second;
630 values[(*iter2).first]+=(*iter2).second;
632 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
633 deno[idx][(*iter2).first]=sum;
636 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
638 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
639 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
643 void MEDCouplingRemapper::computeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
644 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
646 std::map<int,double> values;
648 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
651 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
653 sum+=(*iter2).second;
654 values[(*iter2).first]+=(*iter2).second;
656 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
657 denoReverse[(*iter2).first][idx]=sum;
660 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
662 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
663 deno[idx][(*iter2).first]=values[(*iter2).first];
667 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
668 const std::vector< std::map<int,double> >& m2D,
669 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
670 const int *corrCellIdTrg)
672 int nbOf2DCellsTrg=m2D.size();
673 int nbOf1DCellsTrg=m1D.size();
674 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
675 _matrix.resize(nbOf3DCellsTrg);
677 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
679 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
682 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
684 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
686 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
693 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
696 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
698 std::cout << "Target Cell # " << id << " : ";
699 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
700 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
701 std::cout << std::endl;
705 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const
711 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
712 * If not the behaviour is unpredictable.
713 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than \a maxValAbs this coefficient is
714 * 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.
715 * 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
718 * \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.
719 * \return a positive value that tells the number of coefficients put to 0. The 0 returned value means that the matrix has remained unchanged.
720 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix
722 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs(double maxValAbs) throw(INTERP_KERNEL::Exception)
725 std::vector<std::map<int,double> > matrixNew(_matrix.size());
727 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++,i++)
729 std::map<int,double>& rowNew=matrixNew[i];
730 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
732 if(fabs((*it2).second)>maxValAbs)
733 rowNew[(*it2).first]=(*it2).second;
744 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
745 * If not the behaviour is unpredictable.
746 * This method works on precomputed \a this->_matrix. All coefficients in the matrix is lower than delta multiplied by \a scaleFactor this coefficient is
747 * 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.
748 * delta is the value returned by MEDCouplingRemapper::getMaxValueInCrudeMatrix method.
749 * 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
752 * \param [in] scaleFactor is the scale factor from which coefficients lower than \a scaleFactor times range width of coefficients are set to zero.
753 * \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
754 * that all coefficients are null.
755 * \sa MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrixAbs
757 int MEDCouplingRemapper::nullifiedTinyCoeffInCrudeMatrix(double scaleFactor) throw(INTERP_KERNEL::Exception)
759 double maxVal=getMaxValueInCrudeMatrix();
762 return nullifiedTinyCoeffInCrudeMatrixAbs(scaleFactor*maxVal);
766 * This method is supposed to be called , if needed, right after MEDCouplingRemapper::prepare or MEDCouplingRemapper::prepareEx.
767 * If not the behaviour is unpredictable.
768 * This method returns the maximum of the absolute values of coefficients into the sparse crude matrix.
769 * The returned value is positive.
771 double MEDCouplingRemapper::getMaxValueInCrudeMatrix() const throw(INTERP_KERNEL::Exception)
774 for(std::vector<std::map<int,double> >::const_iterator it1=_matrix.begin();it1!=_matrix.end();it1++)
775 for(std::map<int,double>::const_iterator it2=(*it1).begin();it2!=(*it1).end();it2++)
776 if(fabs((*it2).second)>ret)
777 ret=fabs((*it2).second);