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