Back to index

salome-med  6.5.0
MEDMEM_Remapper.cxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #include "MEDMEM_Remapper.hxx"
00021 #include "MEDMEM_Exception.hxx"
00022 #include "Interpolation.hxx"
00023 #include "Interpolation2D.txx"
00024 #include "Interpolation3D.txx"
00025 #include "Interpolation3DSurf.hxx"
00026 
00027 
00028   
00029 //   int InterpolationOptions::_printLevel=0;
00030 //   IntersectionType  InterpolationOptions::_intersectionType=Triangulation;
00031 //   double  InterpolationOptions::_precision=1e-12;;
00032 //   double  InterpolationOptions::_medianPlane =0.5;
00033 //   bool  InterpolationOptions::_doRotate =true;
00034 //   double  InterpolationOptions::_boundingBoxAdjustment =0.1;
00035 //   int  InterpolationOptions::_orientation =0;
00036 //   SplittingPolicy  InterpolationOptions::_splittingPolicy =GENERAL_48;
00037 
00073 MEDMEM_REMAPPER::MEDMEM_REMAPPER():_matrix(0),_sourceMesh(0), _targetMesh(0), _sourceSupport(0), _targetSupport(0)
00074 {
00075 }
00076 
00077 MEDMEM_REMAPPER::~MEDMEM_REMAPPER()
00078 {
00079   delete _matrix;
00080   if(_sourceMesh)
00081     _sourceMesh->removeReference();
00082   if(_targetMesh)
00083     _targetMesh->removeReference();
00084 //   if(_sourceSupport)
00085 //     _sourceSupport->removeReference();
00086 //   if(_targetSupport)
00087 //     _targetSupport->removeReference();
00088 }
00095 int MEDMEM_REMAPPER::prepare(const MEDMEM::MESH& mesh_source, const MEDMEM::MESH& mesh_target, const char *method)
00096 {
00097   const int sm_spacedim = mesh_source.getSpaceDimension();
00098   const int tm_spacedim = mesh_target.getSpaceDimension();
00099   int sm_meshdim = mesh_source.getMeshDimension();
00100   int tm_meshdim = mesh_target.getMeshDimension();
00101 
00102   delete _matrix;
00103   _matrix= new INTERP_KERNEL::Matrix<double,INTERP_KERNEL::ALL_FORTRAN_MODE>;
00104 
00105   if(_sourceMesh)
00106     _sourceMesh->removeReference();
00107   _sourceMesh= new MEDMEM::MESH((MEDMEM::MESH&)mesh_source);
00108   if(_targetMesh)
00109     _targetMesh->removeReference();
00110   _targetMesh= new MEDMEM::MESH((MEDMEM::MESH&)mesh_target);
00111 
00112   std::string methodC=method;
00113   if(methodC == "P0P0"){
00114     _sourceFieldType = "P0";
00115     _targetFieldType = "P0";
00116   }
00117   else if(methodC =="P0P1"){
00118     _sourceFieldType = "P0";
00119     _targetFieldType = "P1";
00120   }
00121   else if(methodC == "P1P0"){
00122     _sourceFieldType = "P1";
00123     _targetFieldType = "P0";
00124   }
00125   else if(methodC == "P1P1"){
00126     _sourceFieldType = "P1";
00127     _targetFieldType = "P1";
00128   }
00129   else
00130     throw INTERP_KERNEL::Exception("MEDMEM_REMAPPER::prepare: Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
00131                 
00132 
00133 //   if(_sourceSupport)
00134 //     _sourceSupport->removeReference();
00135 //   if(_targetSupport)
00136 //     _targetSupport->removeReference();
00137   if(   _sourceFieldType == "P0")
00138     _sourceSupport = ((MEDMEM::MESH *)_sourceMesh)->getSupportOnAll(MED_EN::MED_CELL);
00139   else
00140     _sourceSupport = ((MEDMEM::MESH *)_sourceMesh)->getSupportOnAll(MED_EN::MED_NODE);
00141   if(   _targetFieldType == "P0")
00142     _targetSupport = ((MEDMEM::MESH *)_targetMesh)->getSupportOnAll(MED_EN::MED_CELL);
00143   else
00144     _targetSupport = ((MEDMEM::MESH *)_targetMesh)->getSupportOnAll(MED_EN::MED_NODE);
00145         
00146   if (tm_spacedim!=sm_spacedim || tm_meshdim!=sm_meshdim)
00147     throw MEDEXCEPTION("incompatible mesh and/or space dimensions in meshes");
00148   if ((sm_spacedim==2)&&(sm_meshdim==2))
00149     {
00150       MEDNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(_sourceMesh);
00151       MEDNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(_targetMesh);
00152       INTERP_KERNEL::Interpolation2D interpolation(*this);
00153       _nb_cols = interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
00154     }
00155   else if ((sm_spacedim==3)&&(sm_meshdim==3))
00156     {
00157       MEDNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(_sourceMesh);
00158       MEDNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(_targetMesh);
00159       INTERP_KERNEL::Interpolation3D interpolation(*this);
00160       _nb_cols = interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
00161     }
00162   else if ((sm_spacedim==3)&&(sm_meshdim==2))
00163     {
00164       MEDNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(_sourceMesh);
00165       MEDNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(_targetMesh);
00166       INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
00167       _nb_cols =  interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
00168     }
00169   else
00170     throw MEDEXCEPTION("no Interpolation exists for the given mesh and space dimensions");
00171 
00172   _nb_rows=(*_matrix).getNbRows();
00173 
00174   _deno_multiply.resize(_nb_rows);
00175   _matrix->rowSum(_deno_multiply);
00176   _deno_reverse_multiply.resize(_nb_cols);
00177   _matrix->colSum(_deno_reverse_multiply, _nb_cols);
00178 
00179   return 1;
00180 }
00181 
00182 
00183 /*
00184   remaps a vector source field defined on the source mesh onto the target mesh using the intersection matrix
00185   \param field_source : source field to be remapped
00186   \param target_field : field resulting from the remapping on the target mesh
00187 */
00188 void MEDMEM_REMAPPER::transfer(const MEDMEM::FIELD<double>& field_source, MEDMEM::FIELD<double>& field_target)
00189 {
00190   int source_nbcomp=field_source.getNumberOfComponents();
00191   int nb_source_values= field_source.getNumberOfValues();
00192   const double* value_source = field_source.getValue();
00193 
00194   int target_nbcomp=field_target.getNumberOfComponents();
00195   double* value_target = const_cast<double*> (field_target.getValue());
00196 
00197   double precision = getPrecision();
00198 
00199   if (source_nbcomp != target_nbcomp)
00200     throw MEDMEM::MEDEXCEPTION("MEDMEM_REMAPPER::transfer: incoherent number of components for source and target fields");
00201   if (_nb_cols != nb_source_values)
00202     throw MEDMEM::MEDEXCEPTION("MEDMEM_REMAPPER::transfer: incoherent number of field values, cannot cannot multiply by interpolation matrix");
00203 
00204   _matrix->multiply(value_source, value_target,source_nbcomp);
00205 
00206   for (int i=0; i< _nb_rows; i++)
00207     if(fabs(_deno_multiply[i]) > precision)
00208       for(int comp = 0; comp < source_nbcomp; comp++)
00209         value_target[i*source_nbcomp+comp]/=_deno_multiply[i];
00210     else
00211       for(int comp = 0; comp < source_nbcomp; comp++)
00212         value_target[i*source_nbcomp+comp]=0.;            
00213 }
00214 
00215 /*
00216   reverses the direct transfer remapping: a Vector field supported on the target mesh is remapped onto 
00217   the source mesh using the transpose of the intersection matrix
00218   \param field_target : target field to be remapped
00219   \param source_field : field resulting from the remapping on the source mesh
00220 */
00221 void MEDMEM_REMAPPER::reverseTransfer( MEDMEM::FIELD<double>& field_source, const MEDMEM::FIELD<double>& field_target)
00222 {
00223   int source_nbcomp=field_source.getNumberOfComponents();
00224   double* value_source = const_cast<double*> (field_source.getValue());
00225 
00226   int target_nbcomp=field_target.getNumberOfComponents();
00227   int nb_target_values= field_target.getNumberOfValues();
00228   const double* value_target = field_target.getValue();
00229 
00230   double precision = getPrecision();
00231 
00232   if (source_nbcomp != target_nbcomp)
00233     throw MEDMEM::MEDEXCEPTION(" MEDMEM_REMAPPER::reverseTransfer: incoherent number of components for source and target fields");
00234   if (_nb_rows != nb_target_values)
00235     throw MEDMEM::MEDEXCEPTION(" MEDMEM_REMAPPER::reverseTransfer: incoherent number of field values, cannot cannot multiply by interpolation matrix");
00236 
00237   _matrix->transposeMultiply(value_target, value_source, _nb_cols,target_nbcomp);//transposeMultiply(input,output, nbcols,nbcomp)
00238 
00239   for (int i=0; i< _nb_cols; i++)
00240     if(fabs(_deno_reverse_multiply[i]) > precision)
00241       for(int comp = 0; comp < target_nbcomp; comp++)
00242         value_source[i*target_nbcomp+comp]/=_deno_reverse_multiply[i];
00243     else
00244       for(int comp = 0; comp < target_nbcomp; comp++)
00245         value_source[i*target_nbcomp+comp]=0.;
00246 }
00247 
00248 /*
00249   remaps a vector source field defined on the source mesh onto the target mesh using the intersection matrix
00250   \param field_source : source field to be remapped
00251   \return : field resulting from the remapping on the target mesh
00252 */
00253 MEDMEM::FIELD<double> *  MEDMEM_REMAPPER::transferField(const MEDMEM::FIELD<double>& field_source)
00254 {
00255   int source_nbcomp=field_source.getNumberOfComponents();
00256   const double* value_source = field_source.getValue();
00257   int nb_source_values= field_source.getNumberOfValues();
00258 
00259   double precision = getPrecision();
00260 
00261   if (_nb_cols != nb_source_values)
00262     throw MEDMEM::MEDEXCEPTION("MEDMEM_REMAPPER::transfer: incoherent number of field values, cannot cannot multiply by interpolation matrix");
00263         
00264   MEDMEM::FIELD<double> * target_field = new MEDMEM::FIELD<double>(_targetSupport,source_nbcomp);
00265   double* value_target = const_cast<double*> ((*target_field).getValue());
00266                 
00267   _matrix->multiply(value_source, value_target,source_nbcomp);
00268                 
00269   for (int i=0; i< _nb_rows; i++)
00270     if(fabs(_deno_multiply[i]) > precision)
00271       for(int comp = 0; comp < source_nbcomp; comp++)
00272         value_target[i*source_nbcomp+comp]/=_deno_multiply[i];    
00273     else
00274       for(int comp = 0; comp < source_nbcomp; comp++)
00275         value_target[i*source_nbcomp+comp]=0.;
00276         
00277   return target_field;
00278 }
00279 
00280 /*
00281   reverses the direct transfer remapping: a Vector field supported on the target mesh is remapped onto 
00282   the source mesh using the transpose of the intersection matrix
00283   \param field_target : target field to be remapped
00284   \return : field resulting from the remapping on the source mesh
00285 */
00286 MEDMEM::FIELD<double> *  MEDMEM_REMAPPER::reverseTransferField(const MEDMEM::FIELD<double>& field_target)
00287 {
00288   int target_nbcomp=field_target.getNumberOfComponents();
00289   const double* value_target = field_target.getValue();
00290   int nb_target_values= field_target.getNumberOfValues();
00291 
00292   double precision = getPrecision();
00293 
00294   if (_nb_rows != nb_target_values)
00295     throw MEDMEM::MEDEXCEPTION(" MEDMEM_REMAPPER::reverseTransfer: incoherent number of field values, cannot cannot transpose-multiply by interpolation matrix");
00296         
00297   MEDMEM::FIELD<double> * source_field = new MEDMEM::FIELD<double>(_sourceSupport,target_nbcomp);
00298   double* value_source = const_cast<double*> ((*source_field).getValue());
00299         
00300   _matrix->transposeMultiply(value_target, value_source, _nb_cols,target_nbcomp);//transposeMultiply(input,output, nbcols,nbcomp)
00301         
00302   for (int i=0; i< _nb_cols; i++)
00303     if(fabs(_deno_reverse_multiply[i]) > precision)
00304       for(int comp = 0; comp < target_nbcomp; comp++)
00305         value_source[i*target_nbcomp+comp]/=_deno_reverse_multiply[i];
00306     else
00307       for(int comp = 0; comp < target_nbcomp; comp++)
00308         value_source[i*target_nbcomp+comp]=0.;
00309         
00310   return source_field;
00311 }
00312 
00313 int MEDMEM_REMAPPER::setOptionDouble(const std::string& key, double value)
00314 {
00315   bool result = INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
00316         
00317   if(result)
00318     return 1;
00319   else
00320     return 0;
00321 }
00322 
00323 int MEDMEM_REMAPPER::setOptionInt(const std::string& key, int value)
00324 {
00325   bool result = INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
00326         
00327   if(result)
00328     return 1;
00329   else
00330     return 0;
00331 }
00332 
00333 int MEDMEM_REMAPPER::setOptionString(const std::string& key, std::string& value)
00334 {
00335   bool result = INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
00336         
00337   if(result)
00338     return 1;
00339   else
00340     return 0;
00341 }
00342 
00352 MEDMEM::FIELD<double>* MEDMEM_REMAPPER::getSupportVolumes(const MEDMEM::SUPPORT& support)
00353 {
00354   const MEDMEM::GMESH* mesh=support.getMesh();
00355   int dim = mesh->getMeshDimension();
00356   switch (dim)
00357     {
00358     case 2:
00359       return mesh->getArea(&support);
00360     case 3:
00361       return mesh->getVolume(&support);
00362     default:
00363       throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension");
00364     }
00365 }
00366 
00367 void MEDMEM_REMAPPER::printMatrixInfo()
00368 {
00369         std::cout << *_matrix << std::endl;
00370 }
00371