Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2016-2018 The VES code team
       3             :    (see the PEOPLE-VES file at the root of this folder for a list of names)
       4             : 
       5             :    See http://www.ves-code.org for more information.
       6             : 
       7             :    This file is part of VES code module.
       8             : 
       9             :    The VES code module is free software: you can redistribute it and/or modify
      10             :    it under the terms of the GNU Lesser General Public License as published by
      11             :    the Free Software Foundation, either version 3 of the License, or
      12             :    (at your option) any later version.
      13             : 
      14             :    The VES code module is distributed in the hope that it will be useful,
      15             :    but WITHOUT ANY WARRANTY; without even the implied warranty of
      16             :    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
      17             :    GNU Lesser General Public License for more details.
      18             : 
      19             :    You should have received a copy of the GNU Lesser General Public License
      20             :    along with the VES code module.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : 
      23             : #include "GridLinearInterpolation.h"
      24             : 
      25             : #include "tools/Grid.h"
      26             : #include "tools/Exception.h"
      27             : #include "tools/Tools.h"
      28             : 
      29             : 
      30             : namespace PLMD {
      31             : namespace ves {
      32             : 
      33             : 
      34        1869 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_1D(Grid* grid_pntr, const std::vector<double>& arg) {
      35             : 
      36        1869 :   plumed_massert(grid_pntr->getDimension()==1,"The grid is of the wrong dimension, should be one-dimensional");
      37        1869 :   plumed_massert(arg.size()==1,"input value is of the wrong size");
      38             : 
      39        1869 :   double x = arg[0];
      40        5607 :   double grid_dx = grid_pntr->getDx()[0];
      41        3738 :   double grid_min; Tools::convert( grid_pntr->getMin()[0], grid_min);
      42        3738 :   std::vector<unsigned int> i0(1); i0[0] = unsigned( std::floor( (x-grid_min)/grid_dx ) );
      43        3738 :   std::vector<unsigned int> i1(1); i1[0] = unsigned( std::ceil(  (x-grid_min)/grid_dx ) );
      44             :   //
      45        5607 :   double x0 = grid_pntr->getPoint(i0)[0];
      46        5607 :   double x1 = grid_pntr->getPoint(i1)[0];
      47        1869 :   double y0 = grid_pntr->getValue(i0);
      48        1869 :   double y1 = grid_pntr->getValue(i1);
      49             :   //
      50        1869 :   return linearInterpolation(x,x0,x1,y0,y1);
      51             : }
      52             : 
      53             : 
      54       47164 : double GridLinearInterpolation::getGridValueWithLinearInterpolation_2D(Grid* grid_pntr, const std::vector<double>& arg) {
      55       47164 :   plumed_massert(grid_pntr->getDimension()==2,"The grid is of the wrong dimension, should be two-dimensional");
      56       47164 :   plumed_massert(arg.size()==2,"input value is of the wrong size");
      57             : 
      58       47164 :   std::vector<double> grid_delta = grid_pntr->getDx();
      59       47164 :   std::vector<double> grid_min(2);
      60       94328 :   Tools::convert( grid_pntr->getMin()[0], grid_min[0]);
      61       94328 :   Tools::convert( grid_pntr->getMin()[1], grid_min[1]);
      62             : 
      63       47164 :   std::vector<unsigned int> i00(2);
      64       47164 :   std::vector<unsigned int> i01(2);
      65       47164 :   std::vector<unsigned int> i10(2);
      66       47164 :   std::vector<unsigned int> i11(2);
      67             : 
      68      235820 :   i00[0] = i01[0] = unsigned( std::floor( (arg[0]-grid_min[0])/grid_delta[0] ) );
      69      235820 :   i10[0] = i11[0] = unsigned( std::ceil(  (arg[0]-grid_min[0])/grid_delta[0] ) );
      70             : 
      71      235820 :   i00[1] = i10[1] = unsigned( std::floor( (arg[1]-grid_min[1])/grid_delta[1] ) );
      72      235820 :   i01[1] = i11[1] = unsigned( std::ceil(  (arg[1]-grid_min[1])/grid_delta[1] ) );
      73             : 
      74             :   // https://en.wikipedia.org/wiki/Bilinear_interpolation
      75       47164 :   double x = arg[0];
      76       47164 :   double y = arg[1];
      77             : 
      78      141492 :   double x0 = grid_pntr->getPoint(i00)[0];
      79      141492 :   double x1 = grid_pntr->getPoint(i10)[0];
      80      141492 :   double y0 = grid_pntr->getPoint(i00)[1];
      81      141492 :   double y1 = grid_pntr->getPoint(i11)[1];
      82             : 
      83       47164 :   double f00 = grid_pntr->getValue(i00);
      84       47164 :   double f01 = grid_pntr->getValue(i01);
      85       47164 :   double f10 = grid_pntr->getValue(i10);
      86       47164 :   double f11 = grid_pntr->getValue(i11);
      87             : 
      88             :   // linear interpolation in x-direction
      89             :   double fx0 = linearInterpolation(x,x0,x1,f00,f10);
      90             :   double fx1 = linearInterpolation(x,x0,x1,f01,f11);
      91             :   // linear interpolation in y-direction
      92             :   double fxy = linearInterpolation(y,y0,y1,fx0,fx1);
      93             :   //
      94       47164 :   return fxy;
      95             : }
      96             : 
      97             : 
      98       49033 : double GridLinearInterpolation::getGridValueWithLinearInterpolation(Grid* grid_pntr, const std::vector<double>& arg) {
      99       49033 :   unsigned int dim = grid_pntr->getDimension();
     100       49033 :   if(dim==1) {
     101        1869 :     return getGridValueWithLinearInterpolation_1D(grid_pntr,arg);
     102             :   }
     103       47164 :   else if(dim==2) {
     104       47164 :     return getGridValueWithLinearInterpolation_2D(grid_pntr,arg);
     105             :   }
     106             :   else {
     107           0 :     plumed_merror("GridLinearInterpolation::getGridValueWithLinearInterpolation only defined for 1D or 2D grids");
     108             :   }
     109             : }
     110             : 
     111             : 
     112             : 
     113             : }
     114             : }