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
20 #include "MEDCouplingRemapper.hxx"
21 #include "MEDCouplingMemArray.hxx"
22 #include "MEDCouplingFieldDouble.hxx"
23 #include "MEDCouplingFieldTemplate.hxx"
24 #include "MEDCouplingFieldDiscretization.hxx"
25 #include "MEDCouplingExtrudedMesh.hxx"
26 #include "MEDCouplingNormalizedUnstructuredMesh.txx"
28 #include "Interpolation1D.txx"
29 #include "Interpolation2DCurve.hxx"
30 #include "Interpolation2D.txx"
31 #include "Interpolation3D.txx"
32 #include "Interpolation3DSurf.hxx"
33 #include "Interpolation2D1D.txx"
34 #include "Interpolation3D2D.txx"
36 using namespace ParaMEDMEM;
38 MEDCouplingRemapper::MEDCouplingRemapper():_src_mesh(0),_target_mesh(0),_nature_of_deno(NoNature),_time_deno_update(0)
42 MEDCouplingRemapper::~MEDCouplingRemapper()
47 int MEDCouplingRemapper::prepare(const MEDCouplingMesh *srcMesh, const MEDCouplingMesh *targetMesh, const char *method) throw(INTERP_KERNEL::Exception)
50 _src_mesh=const_cast<MEDCouplingMesh *>(srcMesh); _target_mesh=const_cast<MEDCouplingMesh *>(targetMesh);
51 _src_mesh->incrRef(); _target_mesh->incrRef();
52 int meshInterpType=((int)_src_mesh->getType()*16)+(int)_target_mesh->getType();
53 switch(meshInterpType)
55 case 85://Unstructured-Unstructured
56 return prepareUU(method);
57 case 136://Extruded-Extruded
58 return prepareEE(method);
60 throw INTERP_KERNEL::Exception("Not managed type of meshes !");
64 int MEDCouplingRemapper::prepareEx(const MEDCouplingFieldTemplate *src, const MEDCouplingFieldTemplate *target) throw(INTERP_KERNEL::Exception)
66 std::string meth(src->getDiscretization()->getStringRepr());
67 meth+=target->getDiscretization()->getStringRepr();
68 return prepare(src->getMesh(),target->getMesh(),meth.c_str());
72 * This method performs the operation source to target using matrix computed in ParaMEDMEM::MEDCouplingRemapper::prepare method.
73 * 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.
75 * \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.
76 * \param [out] targetField the destination field with the allocated array in which all tuples will be overwritten.
78 void MEDCouplingRemapper::transfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
80 transferUnderground(srcField,targetField,true,dftValue);
84 * This method is equivalent to ParaMEDMEM::MEDCouplingRemapper::transfer except that here \b targetField is a in/out parameter.
85 * 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
87 * This method requires that \b targetField was fully defined and allocated. If the array is not allocated an exception will be thrown.
89 * \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.
90 * \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.
92 void MEDCouplingRemapper::partialTransfer(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField) throw(INTERP_KERNEL::Exception)
94 transferUnderground(srcField,targetField,false,std::numeric_limits<double>::max());
97 void MEDCouplingRemapper::reverseTransfer(MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
99 if(_src_method!=srcField->getDiscretization()->getStringRepr())
100 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
101 if(_target_method!=targetField->getDiscretization()->getStringRepr())
102 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
103 if(srcField->getNature()!=targetField->getNature())
104 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
105 DataArrayDouble *array=srcField->getArray();
106 int trgNbOfCompo=targetField->getNumberOfComponents();
109 if(trgNbOfCompo!=srcField->getNumberOfComponents())
110 throw INTERP_KERNEL::Exception("Number of components mismatch !");
114 array=DataArrayDouble::New();
115 array->alloc(srcField->getNumberOfTuples(),trgNbOfCompo);
116 srcField->setArray(array);
119 computeDeno(srcField->getNature(),srcField,targetField);
120 double *resPointer=array->getPointer();
121 const double *inputPointer=targetField->getArray()->getConstPointer();
122 computeReverseProduct(inputPointer,trgNbOfCompo,dftValue,resPointer);
125 MEDCouplingFieldDouble *MEDCouplingRemapper::transferField(const MEDCouplingFieldDouble *srcField, double dftValue) throw(INTERP_KERNEL::Exception)
127 if(_src_method!=srcField->getDiscretization()->getStringRepr())
128 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
129 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),srcField->getTimeDiscretization());
130 ret->copyTinyAttrFrom(srcField);
131 ret->setNature(srcField->getNature());
132 ret->setMesh(_target_mesh);
133 transfer(srcField,ret,dftValue);
137 MEDCouplingFieldDouble *MEDCouplingRemapper::reverseTransferField(const MEDCouplingFieldDouble *targetField, double dftValue) throw(INTERP_KERNEL::Exception)
139 if(_target_method!=targetField->getDiscretization()->getStringRepr())
140 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
141 MEDCouplingFieldDouble *ret=MEDCouplingFieldDouble::New(MEDCouplingFieldDiscretization::getTypeOfFieldFromStringRepr(_target_method.c_str()),targetField->getTimeDiscretization());
142 ret->copyTinyAttrFrom(targetField);
143 ret->setNature(targetField->getNature());
144 ret->setMesh(_src_mesh);
145 reverseTransfer(ret,targetField,dftValue);
149 bool MEDCouplingRemapper::setOptionInt(const std::string& key, int value)
151 return INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
154 bool MEDCouplingRemapper::setOptionDouble(const std::string& key, double value)
156 return INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
159 bool MEDCouplingRemapper::setOptionString(const std::string& key, const std::string& value)
161 return INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
164 int MEDCouplingRemapper::prepareUU(const char *method) throw(INTERP_KERNEL::Exception)
166 MEDCouplingUMesh *src_mesh=(MEDCouplingUMesh *)_src_mesh;
167 MEDCouplingUMesh *target_mesh=(MEDCouplingUMesh *)_target_mesh;
168 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
169 const int srcMeshDim=src_mesh->getMeshDimension();
172 srcSpaceDim=src_mesh->getSpaceDimension();
173 const int trgMeshDim=target_mesh->getMeshDimension();
176 trgSpaceDim=target_mesh->getSpaceDimension();
177 if(trgSpaceDim!=srcSpaceDim)
178 if(trgSpaceDim!=-1 && srcSpaceDim!=-1)
179 throw INTERP_KERNEL::Exception("Incoherent space dimension detected between target and source.");
181 if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==1)
183 MEDCouplingNormalizedUnstructuredMesh<1,1> source_mesh_wrapper(src_mesh);
184 MEDCouplingNormalizedUnstructuredMesh<1,1> target_mesh_wrapper(target_mesh);
185 INTERP_KERNEL::Interpolation1D interpolation(*this);
186 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
188 else if(srcMeshDim==1 && trgMeshDim==1 && srcSpaceDim==2)
190 MEDCouplingNormalizedUnstructuredMesh<2,1> source_mesh_wrapper(src_mesh);
191 MEDCouplingNormalizedUnstructuredMesh<2,1> target_mesh_wrapper(target_mesh);
192 INTERP_KERNEL::Interpolation2DCurve interpolation(*this);
193 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
195 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==2)
197 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
198 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
199 INTERP_KERNEL::Interpolation2D interpolation(*this);
200 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
202 else if(srcMeshDim==3 && trgMeshDim==3 && srcSpaceDim==3)
204 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
205 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
206 INTERP_KERNEL::Interpolation3D interpolation(*this);
207 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
209 else if(srcMeshDim==2 && trgMeshDim==2 && srcSpaceDim==3)
211 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
212 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh);
213 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
214 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
216 else if(srcMeshDim==3 && trgMeshDim==1 && srcSpaceDim==3)
218 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
219 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
220 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
221 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
222 INTERP_KERNEL::Interpolation3D interpolation(*this);
223 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
225 else if(srcMeshDim==1 && trgMeshDim==3 && srcSpaceDim==3)
227 if(getIntersectionType()!=INTERP_KERNEL::PointLocator)
228 throw INTERP_KERNEL::Exception("Invalid interpolation requested between 3D and 1D ! Select PointLocator as intersection type !");
229 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
230 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
231 INTERP_KERNEL::Interpolation3D interpolation(*this);
232 std::vector<std::map<int,double> > matrixTmp;
233 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
234 reverseMatrix(matrixTmp,nbCols,_matrix);
235 nbCols=matrixTmp.size();
237 else if(srcMeshDim==2 && trgMeshDim==1 && srcSpaceDim==2)
239 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
241 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
242 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
243 INTERP_KERNEL::Interpolation2D interpolation(*this);
244 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
248 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
249 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
250 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
251 std::vector<std::map<int,double> > matrixTmp;
252 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
253 reverseMatrix(matrixTmp,nbCols,_matrix);
254 nbCols=matrixTmp.size();
255 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
256 if(!duplicateFaces.empty())
258 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
259 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
261 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
262 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
268 else if(srcMeshDim==1 && trgMeshDim==2 && srcSpaceDim==2)
270 if(getIntersectionType()==INTERP_KERNEL::PointLocator)
272 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
273 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
274 INTERP_KERNEL::Interpolation2D interpolation(*this);
275 std::vector<std::map<int,double> > matrixTmp;
276 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
277 reverseMatrix(matrixTmp,nbCols,_matrix);
278 nbCols=matrixTmp.size();
282 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
283 MEDCouplingNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(target_mesh);
284 INTERP_KERNEL::Interpolation2D1D interpolation(*this);
285 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
286 INTERP_KERNEL::Interpolation2D1D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
287 if(!duplicateFaces.empty())
289 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 1D Cells are part of edges shared by 2D cells :\n";
290 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
292 oss << "1D Cell #" << (*it).first << " is part of common edge of following 2D cells ids : ";
293 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
299 else if(srcMeshDim==2 && trgMeshDim==3 && srcSpaceDim==3)
301 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
302 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
303 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
304 nbCols=interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,_matrix,method);
305 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
306 if(!duplicateFaces.empty())
308 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
309 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
311 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
312 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
317 else if(srcMeshDim==3 && trgMeshDim==2 && srcSpaceDim==3)
319 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
320 MEDCouplingNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(target_mesh);
321 INTERP_KERNEL::Interpolation3D2D interpolation(*this);
322 std::vector<std::map<int,double> > matrixTmp;
323 nbCols=interpolation.interpolateMeshes(target_mesh_wrapper,source_mesh_wrapper,matrixTmp,method);
324 reverseMatrix(matrixTmp,nbCols,_matrix);
325 nbCols=matrixTmp.size();
326 INTERP_KERNEL::Interpolation3D2D::DuplicateFacesType duplicateFaces=interpolation.retrieveDuplicateFaces();
327 if(!duplicateFaces.empty())
329 std::ostringstream oss; oss << "An unexpected situation happend ! For the following 2D Cells are part of edges shared by 3D cells :\n";
330 for(std::map<int,std::set<int> >::const_iterator it=duplicateFaces.begin();it!=duplicateFaces.end();it++)
332 oss << "2D Cell #" << (*it).first << " is part of common face of following 3D cells ids : ";
333 std::copy((*it).second.begin(),(*it).second.end(),std::ostream_iterator<int>(oss," "));
338 else if(trgMeshDim==-1)
340 if(srcMeshDim==2 && srcSpaceDim==2)
342 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(src_mesh);
343 INTERP_KERNEL::Interpolation2D interpolation(*this);
344 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
346 else if(srcMeshDim==3 && srcSpaceDim==3)
348 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(src_mesh);
349 INTERP_KERNEL::Interpolation3D interpolation(*this);
350 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
352 else if(srcMeshDim==2 && srcSpaceDim==3)
354 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh);
355 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
356 nbCols=interpolation.toIntegralUniform(source_mesh_wrapper,_matrix,_src_method.c_str());
359 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh to -1D targetMesh");
361 else if(srcMeshDim==-1)
363 if(trgMeshDim==2 && trgSpaceDim==2)
365 MEDCouplingNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(target_mesh);
366 INTERP_KERNEL::Interpolation2D interpolation(*this);
367 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
369 else if(trgMeshDim==3 && trgSpaceDim==3)
371 MEDCouplingNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(target_mesh);
372 INTERP_KERNEL::Interpolation3D interpolation(*this);
373 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
375 else if(trgMeshDim==2 && trgSpaceDim==3)
377 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(target_mesh);
378 INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
379 nbCols=interpolation.fromIntegralUniform(source_mesh_wrapper,_matrix,_target_method.c_str());
382 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension of source mesh from -1D sourceMesh");
385 throw INTERP_KERNEL::Exception("No interpolation available for the given mesh and space dimension");
386 _deno_multiply.clear();
387 _deno_multiply.resize(_matrix.size());
388 _deno_reverse_multiply.clear();
389 _deno_reverse_multiply.resize(nbCols);
394 int MEDCouplingRemapper::prepareEE(const char *method) throw(INTERP_KERNEL::Exception)
396 MEDCouplingExtrudedMesh *src_mesh=(MEDCouplingExtrudedMesh *)_src_mesh;
397 MEDCouplingExtrudedMesh *target_mesh=(MEDCouplingExtrudedMesh *)_target_mesh;
398 std::string methC(method);
400 throw INTERP_KERNEL::Exception("Only P0P0 method implemented for Extruded/Extruded meshes !");
401 INTERP_KERNEL::Interpolation<INTERP_KERNEL::Interpolation3D>::checkAndSplitInterpolationMethod(method,_src_method,_target_method);
402 MEDCouplingNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(src_mesh->getMesh2D());
403 MEDCouplingNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(target_mesh->getMesh2D());
404 INTERP_KERNEL::Interpolation3DSurf interpolation2D(*this);
405 std::vector<std::map<int,double> > matrix2D;
406 int nbCols2D=interpolation2D.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,matrix2D,method);
407 MEDCouplingUMesh *s1D,*t1D;
409 MEDCouplingExtrudedMesh::Project1DMeshes(src_mesh->getMesh1D(),target_mesh->getMesh1D(),getPrecision(),s1D,t1D,v);
410 MEDCouplingNormalizedUnstructuredMesh<1,1> s1DWrapper(s1D);
411 MEDCouplingNormalizedUnstructuredMesh<1,1> t1DWrapper(t1D);
412 std::vector<std::map<int,double> > matrix1D;
413 INTERP_KERNEL::Interpolation1D interpolation1D(*this);
414 int nbCols1D=interpolation1D.interpolateMeshes(s1DWrapper,t1DWrapper,matrix1D,method);
417 buildFinalInterpolationMatrixByConvolution(matrix1D,matrix2D,src_mesh->getMesh3DIds()->getConstPointer(),nbCols2D,nbCols1D,
418 target_mesh->getMesh3DIds()->getConstPointer());
420 _deno_multiply.clear();
421 _deno_multiply.resize(_matrix.size());
422 _deno_reverse_multiply.clear();
423 _deno_reverse_multiply.resize(nbCols2D*nbCols1D);
428 void MEDCouplingRemapper::updateTime() const
432 void MEDCouplingRemapper::releaseData(bool matrixSuppression)
435 _src_mesh->decrRef();
437 _target_mesh->decrRef();
440 if(matrixSuppression)
443 _deno_multiply.clear();
444 _deno_reverse_multiply.clear();
448 void MEDCouplingRemapper::transferUnderground(const MEDCouplingFieldDouble *srcField, MEDCouplingFieldDouble *targetField, bool isDftVal, double dftValue) throw(INTERP_KERNEL::Exception)
450 if(_src_method!=srcField->getDiscretization()->getStringRepr())
451 throw INTERP_KERNEL::Exception("Incoherency with prepare call for source field");
452 if(_target_method!=targetField->getDiscretization()->getStringRepr())
453 throw INTERP_KERNEL::Exception("Incoherency with prepare call for target field");
454 if(srcField->getNature()!=targetField->getNature())
455 throw INTERP_KERNEL::Exception("Natures of fields mismatch !");
456 DataArrayDouble *array=targetField->getArray();
457 int srcNbOfCompo=srcField->getNumberOfComponents();
460 if(srcNbOfCompo!=targetField->getNumberOfComponents())
461 throw INTERP_KERNEL::Exception("Number of components mismatch !");
466 throw INTERP_KERNEL::Exception("MEDCouplingRemapper::partialTransfer : This method requires that the array of target field exists ! Allocate it or call MEDCouplingRemapper::transfer instead !");
467 array=DataArrayDouble::New();
468 array->alloc(targetField->getNumberOfTuples(),srcNbOfCompo);
469 targetField->setArray(array);
472 computeDeno(srcField->getNature(),srcField,targetField);
473 double *resPointer=array->getPointer();
474 const double *inputPointer=srcField->getArray()->getConstPointer();
475 computeProduct(inputPointer,srcNbOfCompo,isDftVal,dftValue,resPointer);
478 void MEDCouplingRemapper::computeDeno(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField)
481 return computeDenoFromScratch(nat,srcField,trgField);
482 else if(nat!=_nature_of_deno)
483 return computeDenoFromScratch(nat,srcField,trgField);
484 else if(nat==_nature_of_deno && _time_deno_update!=getTimeOfThis())
485 return computeDenoFromScratch(nat,srcField,trgField);
488 void MEDCouplingRemapper::computeDenoFromScratch(NatureOfField nat, const MEDCouplingFieldDouble *srcField, const MEDCouplingFieldDouble *trgField) throw(INTERP_KERNEL::Exception)
491 _time_deno_update=getTimeOfThis();
492 switch(_nature_of_deno)
494 case ConservativeVolumic:
496 computeRowSumAndColSum(_matrix,_deno_multiply,_deno_reverse_multiply);
501 MEDCouplingFieldDouble *deno=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
502 MEDCouplingFieldDouble *denoR=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
503 const double *denoPtr=deno->getArray()->getConstPointer();
504 const double *denoRPtr=denoR->getArray()->getConstPointer();
505 if(trgField->getMesh()->getMeshDimension()==-1)
507 double *denoRPtr2=denoR->getArray()->getPointer();
508 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
510 if(srcField->getMesh()->getMeshDimension()==-1)
512 double *denoPtr2=deno->getArray()->getPointer();
513 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
516 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
517 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
519 _deno_multiply[idx][(*iter2).first]=denoPtr[(*iter2).first];
520 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[idx];
526 case IntegralGlobConstraint:
528 computeColSumAndRowSum(_matrix,_deno_multiply,_deno_reverse_multiply);
533 MEDCouplingFieldDouble *deno=trgField->getDiscretization()->getMeasureField(trgField->getMesh(),true);
534 MEDCouplingFieldDouble *denoR=srcField->getDiscretization()->getMeasureField(srcField->getMesh(),true);
535 const double *denoPtr=deno->getArray()->getConstPointer();
536 const double *denoRPtr=denoR->getArray()->getConstPointer();
537 if(trgField->getMesh()->getMeshDimension()==-1)
539 double *denoRPtr2=denoR->getArray()->getPointer();
540 denoRPtr2[0]=std::accumulate(denoPtr,denoPtr+deno->getNumberOfTuples(),0.);
542 if(srcField->getMesh()->getMeshDimension()==-1)
544 double *denoPtr2=deno->getArray()->getPointer();
545 denoPtr2[0]=std::accumulate(denoRPtr,denoRPtr+denoR->getNumberOfTuples(),0.);
548 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
549 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
551 _deno_multiply[idx][(*iter2).first]=denoPtr[idx];
552 _deno_reverse_multiply[(*iter2).first][idx]=denoRPtr[(*iter2).first];
559 throw INTERP_KERNEL::Exception("No nature specified ! Select one !");
563 void MEDCouplingRemapper::computeProduct(const double *inputPointer, int inputNbOfCompo, bool isDftVal, double dftValue, double *resPointer)
566 double *tmp=new double[inputNbOfCompo];
567 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
572 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
576 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,0.);
577 std::map<int,double>::const_iterator iter3=_deno_multiply[idx].begin();
578 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++,iter3++)
580 std::transform(inputPointer+(*iter2).first*inputNbOfCompo,inputPointer+((*iter2).first+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/(*iter3).second));
581 std::transform(tmp,tmp+inputNbOfCompo,resPointer+idx*inputNbOfCompo,resPointer+idx*inputNbOfCompo,std::plus<double>());
587 void MEDCouplingRemapper::computeReverseProduct(const double *inputPointer, int inputNbOfCompo, double dftValue, double *resPointer)
589 std::vector<bool> isReached(_deno_reverse_multiply.size(),false);
591 double *tmp=new double[inputNbOfCompo];
592 std::fill(resPointer,resPointer+inputNbOfCompo*_deno_reverse_multiply.size(),0.);
593 for(std::vector<std::map<int,double> >::const_iterator iter1=_matrix.begin();iter1!=_matrix.end();iter1++,idx++)
595 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
597 isReached[(*iter2).first]=true;
598 std::transform(inputPointer+idx*inputNbOfCompo,inputPointer+(idx+1)*inputNbOfCompo,tmp,std::bind2nd(std::multiplies<double>(),(*iter2).second/_deno_reverse_multiply[(*iter2).first][idx]));
599 std::transform(tmp,tmp+inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,resPointer+((*iter2).first)*inputNbOfCompo,std::plus<double>());
604 for(std::vector<bool>::const_iterator iter3=isReached.begin();iter3!=isReached.end();iter3++,idx++)
606 std::fill(resPointer+idx*inputNbOfCompo,resPointer+(idx+1)*inputNbOfCompo,dftValue);
609 void MEDCouplingRemapper::reverseMatrix(const std::vector<std::map<int,double> >& matIn, int nbColsMatIn, std::vector<std::map<int,double> >& matOut)
611 matOut.resize(nbColsMatIn);
613 for(std::vector<std::map<int,double> >::const_iterator iter1=matIn.begin();iter1!=matIn.end();iter1++,id++)
614 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
615 matOut[(*iter2).first][id]=(*iter2).second;
618 void MEDCouplingRemapper::computeRowSumAndColSum(const std::vector<std::map<int,double> >& matrixDeno,
619 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
621 std::map<int,double> values;
623 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
626 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
628 sum+=(*iter2).second;
629 values[(*iter2).first]+=(*iter2).second;
631 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
632 deno[idx][(*iter2).first]=sum;
635 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
637 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
638 denoReverse[(*iter2).first][idx]=values[(*iter2).first];
642 void MEDCouplingRemapper::computeColSumAndRowSum(const std::vector<std::map<int,double> >& matrixDeno,
643 std::vector<std::map<int,double> >& deno, std::vector<std::map<int,double> >& denoReverse)
645 std::map<int,double> values;
647 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
650 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
652 sum+=(*iter2).second;
653 values[(*iter2).first]+=(*iter2).second;
655 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
656 denoReverse[(*iter2).first][idx]=sum;
659 for(std::vector<std::map<int,double> >::const_iterator iter1=matrixDeno.begin();iter1!=matrixDeno.end();iter1++,idx++)
661 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
662 deno[idx][(*iter2).first]=values[(*iter2).first];
666 void MEDCouplingRemapper::buildFinalInterpolationMatrixByConvolution(const std::vector< std::map<int,double> >& m1D,
667 const std::vector< std::map<int,double> >& m2D,
668 const int *corrCellIdSrc, int nbOf2DCellsSrc, int nbOf1DCellsSrc,
669 const int *corrCellIdTrg)
671 int nbOf2DCellsTrg=m2D.size();
672 int nbOf1DCellsTrg=m1D.size();
673 int nbOf3DCellsTrg=nbOf2DCellsTrg*nbOf1DCellsTrg;
674 _matrix.resize(nbOf3DCellsTrg);
676 for(std::vector< std::map<int,double> >::const_iterator iter2R=m2D.begin();iter2R!=m2D.end();iter2R++,id2R++)
678 for(std::map<int,double>::const_iterator iter2C=(*iter2R).begin();iter2C!=(*iter2R).end();iter2C++)
681 for(std::vector< std::map<int,double> >::const_iterator iter1R=m1D.begin();iter1R!=m1D.end();iter1R++,id1R++)
683 for(std::map<int,double>::const_iterator iter1C=(*iter1R).begin();iter1C!=(*iter1R).end();iter1C++)
685 _matrix[corrCellIdTrg[id1R*nbOf2DCellsTrg+id2R]][corrCellIdSrc[(*iter1C).first*nbOf2DCellsSrc+(*iter2C).first]]=(*iter1C).second*((*iter2C).second);
692 void MEDCouplingRemapper::PrintMatrix(const std::vector<std::map<int,double> >& m)
695 for(std::vector<std::map<int,double> >::const_iterator iter1=m.begin();iter1!=m.end();iter1++,id++)
697 std::cout << "Target Cell # " << id << " : ";
698 for(std::map<int,double>::const_iterator iter2=(*iter1).begin();iter2!=(*iter1).end();iter2++)
699 std::cout << "(" << (*iter2).first << "," << (*iter2).second << "), ";
700 std::cout << std::endl;
704 const std::vector<std::map<int,double> >& MEDCouplingRemapper::getCrudeMatrix() const