src/MEDMEM/MEDMEM_Remapper.cxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   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.
   7 //
   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.
  12 //
  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
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19
  20 #include "MEDMEM_Remapper.hxx"
  21 #include "MEDMEM_Exception.hxx"
  22 #include "Interpolation.hxx"
  23 #include "Interpolation2D.txx"
  24 #include "Interpolation3D.txx"
  25 #include "Interpolation3DSurf.hxx"
  26
  27
  28
  29 //   int InterpolationOptions::_printLevel=0;
  30 //   IntersectionType  InterpolationOptions::_intersectionType=Triangulation;
  31 //   double  InterpolationOptions::_precision=1e-12;;
  32 //   double  InterpolationOptions::_medianPlane =0.5;
  33 //   bool  InterpolationOptions::_doRotate =true;
  34 //   double  InterpolationOptions::_boundingBoxAdjustment =0.1;
  35 //   int  InterpolationOptions::_orientation =0;
  36 //   SplittingPolicy  InterpolationOptions::_splittingPolicy =GENERAL_48;
  37
  38 /*!
  39     \defgroup medmemremapper MEDMEM::MEDMEM_REMAPPER
  40
  41     \section overview Overview
  42
  43     \c MEDMEM_REMAPPER enables conservative remapping of fields
  44     between two sequential codes.
  45     The computation is possible for 3D meshes and 2D meshes or 3D surfaces. It enables fast sequential localization, based on a bounding box tree structure. It is based on cell-cell intersection or on a point in element search, depending on the field type. Fields can either lie on cells (P0) or nodes (P1).
  46
  47     A typical use of \c MEDMEM_REMAPPER encompasses two distinct phases :
  48     - A setup phase during which the intersection volumes are computed and the interpolation matrix is stored. This corresponds to calling the \c MEDMEM_REMAPPER::prepare() method.
  49     - A use phase during which the remappings are actually performed. This corresponds to the calls to transfer() and reverseTransfer() which actually performs the interpolation.
  50
  51     The following code excerpt illutrates a typical use of the \c MEDMEM_REMAPPER class.
  52
  53     \code
  54     ...
  55     std::string sourceFileName("source.med");
  56     MEDMEM::MESH med_source_mesh(MED_DRIVER,sourceFileName,"Source_Mesh");
  57     std::string targetFileName("target.med");
  58     MEDMEM::MESH med_target_mesh(MED_DRIVER,targetFileName,"Target_Mesh");
  59     FIELD<double> sourceField(MED_DRIVER,sourceFileName,"Density",0,0);
  60     FIELD<double> targetField;
  61     MEDMEM_Remapper mapper;
  62     mapper.setOptionsDouble("Precision",1e-7);
  63     mapper.setOptionsString("Intersection_type",Geometric2D);
  64     mapper.prepare(med_source_mesh,med_target_mesh,"P0P1");
  65     mapper.transfer(sourceField,targetField);
  66     //use targetField
  67     ...
  68     \endcode
  69
  70     @{
  71  */
  72
  73 MEDMEM_REMAPPER::MEDMEM_REMAPPER():_matrix(0),_sourceMesh(0), _targetMesh(0), _sourceSupport(0), _targetSupport(0)
  74 {
  75 }
  76
  77 MEDMEM_REMAPPER::~MEDMEM_REMAPPER()
  78 {
  79   delete _matrix;
  80   if(_sourceMesh)
  81     _sourceMesh->removeReference();
  82   if(_targetMesh)
  83     _targetMesh->removeReference();
  84 //   if(_sourceSupport)
  85 //     _sourceSupport->removeReference();
  86 //   if(_targetSupport)
  87 //     _targetSupport->removeReference();
  88 }
  89 /*! This method computes the intersection matrix between
  90  * source \a mesh_source and \a mesh_target. It is a preliminary step
  91  * that is necessary before calling the \c transfer() method.
  92  * The method analyses the dimensions of the meshes and checks for compatibility.
  93  *
  94  */
  95 int MEDMEM_REMAPPER::prepare(const MEDMEM::MESH& mesh_source, const MEDMEM::MESH& mesh_target, const char *method)
  96 {
  97   const int sm_spacedim = mesh_source.getSpaceDimension();
  98   const int tm_spacedim = mesh_target.getSpaceDimension();
  99   int sm_meshdim = mesh_source.getMeshDimension();
 100   int tm_meshdim = mesh_target.getMeshDimension();
 101
 102   delete _matrix;
 103   _matrix= new INTERP_KERNEL::Matrix<double,INTERP_KERNEL::ALL_FORTRAN_MODE>;
 104
 105   if(_sourceMesh)
 106     _sourceMesh->removeReference();
 107   _sourceMesh= new MEDMEM::MESH((MEDMEM::MESH&)mesh_source);
 108   if(_targetMesh)
 109     _targetMesh->removeReference();
 110   _targetMesh= new MEDMEM::MESH((MEDMEM::MESH&)mesh_target);
 111
 112   std::string methodC=method;
 113   if(methodC == "P0P0"){
 114     _sourceFieldType = "P0";
 115     _targetFieldType = "P0";
 116   }
 117   else if(methodC =="P0P1"){
 118     _sourceFieldType = "P0";
 119     _targetFieldType = "P1";
 120   }
 121   else if(methodC == "P1P0"){
 122     _sourceFieldType = "P1";
 123     _targetFieldType = "P0";
 124   }
 125   else if(methodC == "P1P1"){
 126     _sourceFieldType = "P1";
 127     _targetFieldType = "P1";
 128   }
 129   else
 130     throw INTERP_KERNEL::Exception("MEDMEM_REMAPPER::prepare: Invalid method specified ! Must be in : \"P0P0\" \"P0P1\" \"P1P0\" or \"P1P1\"");
 131
 132
 133 //   if(_sourceSupport)
 134 //     _sourceSupport->removeReference();
 135 //   if(_targetSupport)
 136 //     _targetSupport->removeReference();
 137   if(   _sourceFieldType == "P0")
 138     _sourceSupport = ((MEDMEM::MESH *)_sourceMesh)->getSupportOnAll(MED_EN::MED_CELL);
 139   else
 140     _sourceSupport = ((MEDMEM::MESH *)_sourceMesh)->getSupportOnAll(MED_EN::MED_NODE);
 141   if(   _targetFieldType == "P0")
 142     _targetSupport = ((MEDMEM::MESH *)_targetMesh)->getSupportOnAll(MED_EN::MED_CELL);
 143   else
 144     _targetSupport = ((MEDMEM::MESH *)_targetMesh)->getSupportOnAll(MED_EN::MED_NODE);
 145
 146   if (tm_spacedim!=sm_spacedim || tm_meshdim!=sm_meshdim)
 147     throw MEDEXCEPTION("incompatible mesh and/or space dimensions in meshes");
 148   if ((sm_spacedim==2)&&(sm_meshdim==2))
 149     {
 150       MEDNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(_sourceMesh);
 151       MEDNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(_targetMesh);
 152       INTERP_KERNEL::Interpolation2D interpolation(*this);
 153       _nb_cols = interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
 154     }
 155   else if ((sm_spacedim==3)&&(sm_meshdim==3))
 156     {
 157       MEDNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(_sourceMesh);
 158       MEDNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(_targetMesh);
 159       INTERP_KERNEL::Interpolation3D interpolation(*this);
 160       _nb_cols = interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
 161     }
 162   else if ((sm_spacedim==3)&&(sm_meshdim==2))
 163     {
 164       MEDNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(_sourceMesh);
 165       MEDNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(_targetMesh);
 166       INTERP_KERNEL::Interpolation3DSurf interpolation(*this);
 167       _nb_cols =  interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
 168     }
 169   else
 170     throw MEDEXCEPTION("no Interpolation exists for the given mesh and space dimensions");
 171
 172   _nb_rows=(*_matrix).getNbRows();
 173
 174   _deno_multiply.resize(_nb_rows);
 175   _matrix->rowSum(_deno_multiply);
 176   _deno_reverse_multiply.resize(_nb_cols);
 177   _matrix->colSum(_deno_reverse_multiply, _nb_cols);
 178
 179   return 1;
 180 }
 181
 182
 183 /*
 184   remaps a vector source field defined on the source mesh onto the target mesh using the intersection matrix
 185   \param field_source : source field to be remapped
 186   \param target_field : field resulting from the remapping on the target mesh
 187 */
 188 void MEDMEM_REMAPPER::transfer(const MEDMEM::FIELD<double>& field_source, MEDMEM::FIELD<double>& field_target)
 189 {
 190   int source_nbcomp=field_source.getNumberOfComponents();
 191   int nb_source_values= field_source.getNumberOfValues();
 192   const double* value_source = field_source.getValue();
 193
 194   int target_nbcomp=field_target.getNumberOfComponents();
 195   double* value_target = const_cast<double*> (field_target.getValue());
 196
 197   double precision = getPrecision();
 198
 199   if (source_nbcomp != target_nbcomp)
 200     throw MEDMEM::MEDEXCEPTION("MEDMEM_REMAPPER::transfer: incoherent number of components for source and target fields");
 201   if (_nb_cols != nb_source_values)
 202     throw MEDMEM::MEDEXCEPTION("MEDMEM_REMAPPER::transfer: incoherent number of field values, cannot cannot multiply by interpolation matrix");
 203
 204   _matrix->multiply(value_source, value_target,source_nbcomp);
 205
 206   for (int i=0; i< _nb_rows; i++)
 207     if(fabs(_deno_multiply[i]) > precision)
 208       for(int comp = 0; comp < source_nbcomp; comp++)
 209         value_target[i*source_nbcomp+comp]/=_deno_multiply[i];
 210     else
 211       for(int comp = 0; comp < source_nbcomp; comp++)
 212         value_target[i*source_nbcomp+comp]=0.;
 213 }
 214
 215 /*
 216   reverses the direct transfer remapping: a Vector field supported on the target mesh is remapped onto
 217   the source mesh using the transpose of the intersection matrix
 218   \param field_target : target field to be remapped
 219   \param source_field : field resulting from the remapping on the source mesh
 220 */
 221 void MEDMEM_REMAPPER::reverseTransfer( MEDMEM::FIELD<double>& field_source, const MEDMEM::FIELD<double>& field_target)
 222 {
 223   int source_nbcomp=field_source.getNumberOfComponents();
 224   double* value_source = const_cast<double*> (field_source.getValue());
 225
 226   int target_nbcomp=field_target.getNumberOfComponents();
 227   int nb_target_values= field_target.getNumberOfValues();
 228   const double* value_target = field_target.getValue();
 229
 230   double precision = getPrecision();
 231
 232   if (source_nbcomp != target_nbcomp)
 233     throw MEDMEM::MEDEXCEPTION(" MEDMEM_REMAPPER::reverseTransfer: incoherent number of components for source and target fields");
 234   if (_nb_rows != nb_target_values)
 235     throw MEDMEM::MEDEXCEPTION(" MEDMEM_REMAPPER::reverseTransfer: incoherent number of field values, cannot cannot multiply by interpolation matrix");
 236
 237   _matrix->transposeMultiply(value_target, value_source, _nb_cols,target_nbcomp);//transposeMultiply(input,output, nbcols,nbcomp)
 238
 239   for (int i=0; i< _nb_cols; i++)
 240     if(fabs(_deno_reverse_multiply[i]) > precision)
 241       for(int comp = 0; comp < target_nbcomp; comp++)
 242         value_source[i*target_nbcomp+comp]/=_deno_reverse_multiply[i];
 243     else
 244       for(int comp = 0; comp < target_nbcomp; comp++)
 245         value_source[i*target_nbcomp+comp]=0.;
 246 }
 247
 248 /*
 249   remaps a vector source field defined on the source mesh onto the target mesh using the intersection matrix
 250   \param field_source : source field to be remapped
 251   \return : field resulting from the remapping on the target mesh
 252 */
 253 MEDMEM::FIELD<double> *  MEDMEM_REMAPPER::transferField(const MEDMEM::FIELD<double>& field_source)
 254 {
 255   int source_nbcomp=field_source.getNumberOfComponents();
 256   const double* value_source = field_source.getValue();
 257   int nb_source_values= field_source.getNumberOfValues();
 258
 259   double precision = getPrecision();
 260
 261   if (_nb_cols != nb_source_values)
 262     throw MEDMEM::MEDEXCEPTION("MEDMEM_REMAPPER::transfer: incoherent number of field values, cannot cannot multiply by interpolation matrix");
 263
 264   MEDMEM::FIELD<double> * target_field = new MEDMEM::FIELD<double>(_targetSupport,source_nbcomp);
 265   double* value_target = const_cast<double*> ((*target_field).getValue());
 266
 267   _matrix->multiply(value_source, value_target,source_nbcomp);
 268
 269   for (int i=0; i< _nb_rows; i++)
 270     if(fabs(_deno_multiply[i]) > precision)
 271       for(int comp = 0; comp < source_nbcomp; comp++)
 272         value_target[i*source_nbcomp+comp]/=_deno_multiply[i];
 273     else
 274       for(int comp = 0; comp < source_nbcomp; comp++)
 275         value_target[i*source_nbcomp+comp]=0.;
 276
 277   return target_field;
 278 }
 279
 280 /*
 281   reverses the direct transfer remapping: a Vector field supported on the target mesh is remapped onto
 282   the source mesh using the transpose of the intersection matrix
 283   \param field_target : target field to be remapped
 284   \return : field resulting from the remapping on the source mesh
 285 */
 286 MEDMEM::FIELD<double> *  MEDMEM_REMAPPER::reverseTransferField(const MEDMEM::FIELD<double>& field_target)
 287 {
 288   int target_nbcomp=field_target.getNumberOfComponents();
 289   const double* value_target = field_target.getValue();
 290   int nb_target_values= field_target.getNumberOfValues();
 291
 292   double precision = getPrecision();
 293
 294   if (_nb_rows != nb_target_values)
 295     throw MEDMEM::MEDEXCEPTION(" MEDMEM_REMAPPER::reverseTransfer: incoherent number of field values, cannot cannot transpose-multiply by interpolation matrix");
 296
 297   MEDMEM::FIELD<double> * source_field = new MEDMEM::FIELD<double>(_sourceSupport,target_nbcomp);
 298   double* value_source = const_cast<double*> ((*source_field).getValue());
 299
 300   _matrix->transposeMultiply(value_target, value_source, _nb_cols,target_nbcomp);//transposeMultiply(input,output, nbcols,nbcomp)
 301
 302   for (int i=0; i< _nb_cols; i++)
 303     if(fabs(_deno_reverse_multiply[i]) > precision)
 304       for(int comp = 0; comp < target_nbcomp; comp++)
 305         value_source[i*target_nbcomp+comp]/=_deno_reverse_multiply[i];
 306     else
 307       for(int comp = 0; comp < target_nbcomp; comp++)
 308         value_source[i*target_nbcomp+comp]=0.;
 309
 310   return source_field;
 311 }
 312
 313 int MEDMEM_REMAPPER::setOptionDouble(const std::string& key, double value)
 314 {
 315   bool result = INTERP_KERNEL::InterpolationOptions::setOptionDouble(key,value);
 316
 317   if(result)
 318     return 1;
 319   else
 320     return 0;
 321 }
 322
 323 int MEDMEM_REMAPPER::setOptionInt(const std::string& key, int value)
 324 {
 325   bool result = INTERP_KERNEL::InterpolationOptions::setOptionInt(key,value);
 326
 327   if(result)
 328     return 1;
 329   else
 330     return 0;
 331 }
 332
 333 int MEDMEM_REMAPPER::setOptionString(const std::string& key, std::string& value)
 334 {
 335   bool result = INTERP_KERNEL::InterpolationOptions::setOptionString(key,value);
 336
 337   if(result)
 338     return 1;
 339   else
 340     return 0;
 341 }
 342
 343 /*!
 344   \brief returns the volumes of the cells underlying the support \a support
 345
 346   For 2D geometries, the returned field contains the areas.
 347   For 3D geometries, the returned field contains the volumes.
 348
 349   \param support : support whose cell volumes are required
 350   \return field containing the volumes
 351 */
 352 MEDMEM::FIELD<double>* MEDMEM_REMAPPER::getSupportVolumes(const MEDMEM::SUPPORT& support)
 353 {
 354   const MEDMEM::GMESH* mesh=support.getMesh();
 355   int dim = mesh->getMeshDimension();
 356   switch (dim)
 357     {
 358     case 2:
 359       return mesh->getArea(&support);
 360     case 3:
 361       return mesh->getVolume(&support);
 362     default:
 363       throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension");
 364     }
 365 }
 366
 367 void MEDMEM_REMAPPER::printMatrixInfo()
 368 {
 369         std::cout << *_matrix << std::endl;
 370 }
 371
 372 /*!
 373 @}
 374 */