src/MEDMEM/Remapper.cxx

   1 //  Copyright (C) 2007-2008  CEA/DEN, EDF R&D
   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 #include "Remapper.hxx"
  20 #include "MEDMEM_Exception.hxx"
  21 #include "Interpolation.hxx"
  22 #include "Interpolation2D.txx"
  23 #include "Interpolation3D.txx"
  24 #include "Interpolation3DSurf.txx"
  25 #include "MEDNormalizedUnstructuredMesh.txx"
  26
  27 namespace INTERP_KERNEL
  28 {
  29
  30   //   int InterpolationOptions::_printLevel=0;
  31   //   IntersectionType  InterpolationOptions::_intersectionType=Triangulation;
  32   //   double  InterpolationOptions::_precision=1e-12;;
  33   //   double  InterpolationOptions::_medianPlane =0.5;
  34   //   bool  InterpolationOptions::_doRotate =true;
  35   //   double  InterpolationOptions::_boundingBoxAdjustment =0.1;
  36   //   int  InterpolationOptions::_orientation =0;
  37   //   SplittingPolicy  InterpolationOptions::_splittingPolicy =GENERAL_48;
  38
  39   Remapper::Remapper():_matrix(0)
  40   {
  41   }
  42
  43   Remapper::~Remapper()
  44   {
  45     delete _matrix;
  46   }
  47   /*! This method computes the intersection matrix between
  48    * source \a mesh_source and \a mesh_target. It is a preliminary step
  49    * that is necessary before calling the \c transfer() method.
  50    * The method analyses the dimensions of the meshes and checks for compatibility.
  51    *
  52    */
  53   void Remapper::prepare(const MEDMEM::MESH& mesh_source, const MEDMEM::MESH& mesh_target, const char *method)
  54   {
  55     const int sm_spacedim = mesh_source.getSpaceDimension();
  56     const int tm_spacedim = mesh_target.getSpaceDimension();
  57     int sm_meshdim = mesh_source.getMeshDimension();
  58     int tm_meshdim = mesh_target.getMeshDimension();
  59
  60     if (tm_spacedim!=sm_spacedim || tm_meshdim!=sm_meshdim)
  61       throw MEDEXCEPTION("incompatible mesh and/or space dimensions in meshes");
  62     _matrix= new Matrix<double,ALL_FORTRAN_MODE>;
  63     if ((sm_spacedim==2)&&(sm_meshdim==2))
  64       {
  65         MEDNormalizedUnstructuredMesh<2,2> source_mesh_wrapper(&mesh_source);
  66         MEDNormalizedUnstructuredMesh<2,2> target_mesh_wrapper(&mesh_target);
  67         Interpolation2D interpolation;
  68         interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
  69       }
  70     else if ((sm_spacedim==3)&&(sm_meshdim==3))
  71       {
  72         MEDNormalizedUnstructuredMesh<3,3> source_mesh_wrapper(&mesh_source);
  73         MEDNormalizedUnstructuredMesh<3,3> target_mesh_wrapper(&mesh_target);
  74         Interpolation3D interpolation;
  75         interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
  76       }
  77     else if ((sm_spacedim==3)&&(sm_meshdim==2))
  78       {
  79         MEDNormalizedUnstructuredMesh<3,2> source_mesh_wrapper(&mesh_source);
  80         MEDNormalizedUnstructuredMesh<3,2> target_mesh_wrapper(&mesh_target);
  81         Interpolation3DSurf interpolation;
  82         interpolation.interpolateMeshes(source_mesh_wrapper,target_mesh_wrapper,*_matrix,method);
  83       }
  84     else
  85       throw MEDEXCEPTION("no Interpolation exists for the given mesh and space dimensions");
  86   }
  87
  88   void Remapper::transfer(const MEDMEM::FIELD<double>& field_source,MEDMEM::FIELD<double>& field_target)
  89   {
  90     int source_nbcomp=field_source.getNumberOfComponents();
  91     int target_nbcomp=field_target.getNumberOfComponents();
  92     if (source_nbcomp != target_nbcomp)
  93       throw MEDMEM::MEDEXCEPTION("incoherent number of components for source and target fields");
  94     if (source_nbcomp>1)
  95       throw MEDMEM::MEDEXCEPTION("interpolations with more than one component are not yet handled");
  96     MEDMEM::FIELD<double>* target_volumes = getSupportVolumes(*field_target.getSupport());
  97     int nbelem_target=field_target.getSupport()->getNumberOfElements(MED_EN::MED_ALL_ELEMENTS);
  98
  99     double* value_target = const_cast<double*> (field_target.getValue());
 100     const double* value_source = field_source.getValue();
 101
 102     _matrix->multiply(value_source, value_target);
 103     for (int i=0; i< nbelem_target; i++)
 104       value_target[i]/=target_volumes->getValueIJ(i+1,1);
 105
 106     delete target_volumes;
 107   }
 108
 109   void Remapper::setOptionDouble(const std::string& key, double value)
 110   {
 111   }
 112
 113   void Remapper::setOptionInt(const std::string& key, int value)
 114   {
 115   }
 116
 117   /*!
 118     \brief returns the volumes of the cells underlying the field \a field
 119
 120     For 2D geometries, the returned field contains the areas.
 121     For 3D geometries, the returned field contains the volumes.
 122
 123     \param field field on which cells the volumes are required
 124     \return field containing the volumes
 125   */
 126   MEDMEM::FIELD<double>* Remapper::getSupportVolumes(const MEDMEM::SUPPORT& support)
 127   {
 128     const MEDMEM::MESH* mesh=support.getMesh();
 129     int dim = mesh->getMeshDimension();
 130     switch (dim)
 131       {
 132       case 2:
 133         return mesh->getArea(&support);
 134       case 3:
 135         return mesh->getVolume(&support);
 136       default:
 137         throw MEDMEM::MEDEXCEPTION("interpolation is not available for this dimension");
 138       }
 139   }
 140
 141
 142 }