Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2023 The plumed team
       3             :    (see the PEOPLE file at the root of the distribution for a list of names)
       4             : 
       5             :    See http://www.plumed.org for more information.
       6             : 
       7             :    This file is part of plumed, version 2.
       8             : 
       9             :    plumed 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             :    plumed 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 plumed.  If not, see <http://www.gnu.org/licenses/>.
      21             : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
      22             : #include "Interpolator.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace gridtools {
      26             : 
      27      520407 : double Interpolator::splineInterpolation( const std::vector<double>& x, std::vector<double>& der ) const {
      28      520407 :   return splineInterpolation( View<const double>(x.data(),x.size()), der );
      29             : }
      30             : 
      31      552517 : double Interpolator::splineInterpolation( const View<const double> x, std::vector<double>& der ) const {
      32             :   plumed_dbg_assert( gridobject.getGridType()=="flat" );
      33      552517 :   unsigned dimension = gridobject.getDimension();
      34             : 
      35             :   double X,X2,X3,value=0;
      36      552517 :   der.assign(der.size(),0.0);
      37      552517 :   std::vector<double> fd(dimension);
      38      552517 :   std::vector<double> C(dimension);
      39      552517 :   std::vector<double> D(dimension);
      40      552517 :   std::vector<double> dder(dimension);
      41             : 
      42      552517 :   std::vector<unsigned> nindices(dimension);
      43      552517 :   std::vector<unsigned> indices(dimension);
      44      552517 :   gridobject.getIndices( x, indices );
      45      552517 :   std::vector<double> xfloor(dimension);
      46      552517 :   gridobject.getGridPointCoordinates( gridobject.getIndex(x), nindices, xfloor );
      47             : 
      48             :   // loop over neighbors
      49             :   std::vector<unsigned> neigh;
      50             :   unsigned n_neigh;
      51      552517 :   gridobject.getSplineNeighbors( gridobject.getIndex(indices), n_neigh, neigh );
      52     2738563 :   for(unsigned int ipoint=0; ipoint<n_neigh; ++ipoint) {
      53     2186046 :     double grid=values->get( neigh[ipoint] );
      54     6628656 :     for(unsigned j=0; j<dimension; ++j) {
      55     4442610 :       dder[j] = values->getGridDerivative( neigh[ipoint], j );
      56             :     }
      57             : 
      58     2186046 :     gridobject.getIndices( neigh[ipoint], nindices );
      59             :     double ff=1.0;
      60     6628656 :     for(unsigned j=0; j<dimension; ++j) {
      61             :       int x0=1;
      62     4442610 :       if(nindices[j]==indices[j]) {
      63             :         x0=0;
      64             :       }
      65     4442610 :       double ddx=gridobject.getGridSpacing()[j];
      66     4442610 :       X=fabs((x[j]-xfloor[j])/ddx-(double)x0);
      67     4442610 :       X2=X*X;
      68     4442610 :       X3=X2*X;
      69             :       double yy;
      70     4442610 :       if(fabs(grid)<0.0000001) {
      71             :         yy=0.0;
      72             :       } else {
      73     4415152 :         yy=-dder[j]/grid;
      74             :       }
      75     6684115 :       C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*ddx;
      76     4442610 :       D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*ddx;
      77     4442610 :       D[j]*=(x0?-1.0:1.0)/ddx;
      78     4442610 :       ff*=C[j];
      79             :     }
      80     6628656 :     for(unsigned j=0; j<dimension; ++j) {
      81     4442610 :       fd[j]=D[j];
      82    13614108 :       for(unsigned i=0; i<dimension; ++i)
      83     9171498 :         if(i!=j) {
      84     4728888 :           fd[j]*=C[i];
      85             :         }
      86             :     }
      87     2186046 :     value+=grid*ff;
      88     6628656 :     for(unsigned j=0; j<dimension; ++j) {
      89     4442610 :       der[j]+=grid*fd[j];
      90             :     }
      91             :   }
      92      552517 :   return value;
      93             : }
      94             : 
      95             : }
      96             : }