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      544167 : double Interpolator::splineInterpolation( const std::vector<double>& x, std::vector<double>& der ) const {
      28             :   plumed_dbg_assert( gridobject.getGridType()=="flat" );
      29      544167 :   unsigned dimension = gridobject.getDimension();
      30             : 
      31             :   double X,X2,X3,value=0;
      32      544167 :   der.assign(der.size(),0.0);
      33      544167 :   std::vector<double> fd(dimension);
      34      544167 :   std::vector<double> C(dimension);
      35      544167 :   std::vector<double> D(dimension);
      36      544167 :   std::vector<double> dder(dimension);
      37             : 
      38      544167 :   std::vector<unsigned> nindices(dimension);
      39      544167 :   std::vector<unsigned> indices(dimension);
      40      544167 :   gridobject.getIndices( x, indices );
      41      544167 :   std::vector<double> xfloor(dimension);
      42      544167 :   gridobject.getGridPointCoordinates( gridobject.getIndex(x), nindices, xfloor );
      43             : 
      44             :   // loop over neighbors
      45             :   std::vector<unsigned> neigh;
      46             :   unsigned n_neigh;
      47      544167 :   gridobject.getSplineNeighbors( gridobject.getIndex(indices), n_neigh, neigh );
      48     2713479 :   for(unsigned int ipoint=0; ipoint<n_neigh; ++ipoint) {
      49     2169312 :     double grid=values->get( neigh[ipoint] );
      50     6595096 :     for(unsigned j=0; j<dimension; ++j) {
      51     4425784 :       dder[j] = values->getGridDerivative( neigh[ipoint], j );
      52             :     }
      53             : 
      54     2169312 :     gridobject.getIndices( neigh[ipoint], nindices );
      55             :     double ff=1.0;
      56     6595096 :     for(unsigned j=0; j<dimension; ++j) {
      57             :       int x0=1;
      58     4425784 :       if(nindices[j]==indices[j]) {
      59             :         x0=0;
      60             :       }
      61     4425784 :       double ddx=gridobject.getGridSpacing()[j];
      62     4425784 :       X=fabs((x[j]-xfloor[j])/ddx-(double)x0);
      63     4425784 :       X2=X*X;
      64     4425784 :       X3=X2*X;
      65             :       double yy;
      66     4425784 :       if(fabs(grid)<0.0000001) {
      67             :         yy=0.0;
      68             :       } else {
      69     4398326 :         yy=-dder[j]/grid;
      70             :       }
      71     6658876 :       C[j]=(1.0-3.0*X2+2.0*X3) - (x0?-1.0:1.0)*yy*(X-2.0*X2+X3)*ddx;
      72     4425784 :       D[j]=( -6.0*X +6.0*X2) - (x0?-1.0:1.0)*yy*(1.0-4.0*X +3.0*X2)*ddx;
      73     4425784 :       D[j]*=(x0?-1.0:1.0)/ddx;
      74     4425784 :       ff*=C[j];
      75             :     }
      76     6595096 :     for(unsigned j=0; j<dimension; ++j) {
      77     4425784 :       fd[j]=D[j];
      78    13580176 :       for(unsigned i=0; i<dimension; ++i)
      79     9154392 :         if(i!=j) {
      80     4728608 :           fd[j]*=C[i];
      81             :         }
      82             :     }
      83     2169312 :     value+=grid*ff;
      84     6595096 :     for(unsigned j=0; j<dimension; ++j) {
      85     4425784 :       der[j]+=grid*fd[j];
      86             :     }
      87             :   }
      88      544167 :   return value;
      89             : }
      90             : 
      91             : }
      92             : }