diff --git a/src/tabular_function.cpp b/src/tabular_function.cpp index 3c4922b92b..f14513a41b 100644 --- a/src/tabular_function.cpp +++ b/src/tabular_function.cpp @@ -23,13 +23,6 @@ TabularFunction::TabularFunction() xs(nullptr), ys(nullptr), ys1(nullptr), ys2(nullptr), ys3(nullptr), ys4(nullptr), ys5(nullptr), ys6(nullptr) {} -TabularFunction::TabularFunction(int n, double x1, double x2) - : TabularFunction() -{ - set_xrange(x1, x2); - reset_size(n); -} - TabularFunction::~TabularFunction() { delete [] xs; @@ -50,42 +43,6 @@ void TabularFunction::set_values(int n, double x1, double x2, double *values) initialize(); } -void TabularFunction::set_values(int n, double x1, double x2, double *values, double epsilon) -{ - int ilo=0,ihi=n-1; - double vlo, vhi; - - // get lo/hi boundaries where value < epsilon to shrink stored table - for (int i = ilo; ((fabs(values[i]) <= epsilon) && (i < n)); ++i) - ilo = i; - for (int i = ihi; ((fabs(values[i]) <= epsilon) && (i >= 0)); --i) - ihi = i; - if (ihi < ilo) ihi = ilo; - - vlo = values[ilo]; - vhi = values[ilo]; - for (int i = ilo; i <= ihi; ++i) { - if (vlo > values[i]) vlo = values[i]; - if (vhi < values[i]) vhi = values[i]; - } - - // do not shrink small tables - if (ihi - ilo < 50) { - ilo = 0; - ihi = n-1; - } - - xmin = x1 + (x2-x1)/(n -1)*ilo; - xmax = xmin + (x2-x1)/(n -1)*(ihi-ilo); - xmaxsq = xmax*xmax; - n = ihi - ilo + 1; - reset_size(n); - for (int i = ilo; i <= ihi; i++) { - ys[i-ilo] = values[i]; - } - initialize(); -} - void TabularFunction::set_xrange(double x1, double x2) { xmin = x1; @@ -120,6 +77,10 @@ void TabularFunction::initialize() { int i; rdx = (xmax - xmin) / (size - 1.0); + vmax = 0.0; + for (int i = 0; i < size; i++) { + if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]); + } for (i = 0; i < size; i++) xs[i] = xmin + i * rdx; rdx = 1.0 / rdx; ys1[0] = ys[1] - ys[0]; @@ -140,140 +101,3 @@ void TabularFunction::initialize() ys6[i] = 3.0 * ys3[i] * rdx; } } - -#if 0 - void set_values(int n, double x1, double x2, double * values, double epsilon) - { - int shrink = 1; - int ilo,ihi; - double vlo,vhi; - ilo = 0; - ihi = n-1; - for (int i = 0; i < n; i++) { - if (fabs(values[i]) <= epsilon) { - ilo = i; - } else { - break; - } - } - for (int i = n-1; i >= 0; i--) { - if (fabs(values[i]) <= epsilon) { - ihi = i; - } else { - break; - } - } - if (ihi < ilo) ihi = ilo; - vlo = values[ilo]; - vhi = values[ilo]; - for (int i = ilo; i <= ihi; i++) { - if (vlo > values[i]) vlo = values[i]; - if (vhi < values[i]) vhi = values[i]; - } - -// shrink (remove near zero points) reduces cutoff radius, and therefore computational cost -// do not shrink when x2 < 1.1 (angular function) or x2 > 20.0 (non-radial function) - if (x2 < 1.1 || x2 > 20.0) { - shrink = 0; - } -// do not shrink when when list is abnormally small - if (ihi - ilo < 50) { - shrink = 0; - } -// shrink if it is a constant - if (vhi - vlo <= epsilon) { -// shrink = 1; - } - - if (shrink == 0) { - ilo = 0; - ihi = n-1; - } - xmin = x1 + (x2-x1)/(n -1)*ilo; - xmax = xmin + (x2-x1)/(n -1)*(ihi-ilo); - xmaxsq = xmax*xmax; - n = ihi - ilo + 1; - resize(n); - for (int i = ilo; i <= ihi; i++) { - ys[i-ilo] = values[i]; - } - initialize(); - } - void value(double x, double &y, int ny, double &y1, int ny1) - { - double ps = (x - xmin) * rdx; - int ks = ps + 0.5; - if (ks > size-1) ks = size-1; - if (ks < 0 ) ks = 0; - ps = ps - ks; - if (ny) y = ((ys3[ks]*ps + ys2[ks])*ps + ys1[ks])*ps + ys[ks]; - if (ny1) y1 = (ys6[ks]*ps + ys5[ks])*ps + ys4[ks]; - } - void print_value() - { - printf("%d %f %f %f \n",size,xmin,xmax,rdx); - printf(" \n"); - for (int i = 0; i < size; i++) { - printf("%f %f \n",xs[i],ys[i]); - } - } - - protected: - - void resize(int n) { - if (n != size) { - size = n; - delete [] xs; - xs = new double[n]; - delete [] ys; - ys = new double[n]; - delete [] ys1; - ys1 = new double[n]; - delete [] ys2; - ys2 = new double[n]; - delete [] ys3; - ys3 = new double[n]; - delete [] ys4; - ys4 = new double[n]; - delete [] ys5; - ys5 = new double[n]; - delete [] ys6; - ys6 = new double[n]; - } - } - void initialize() { - int n = size; - rdx = (xmax-xmin)/(n-1.0); - vmax = 0.0; - for (int i = 0; i < n; i++) { - if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]); - } - for (int i = 0; i < n; i++) { - xs[i] = xmin+i*rdx; - } - rdx = 1.0 / rdx; - ys1[0] = ys[1] - ys[0]; - ys1[1] = 0.5 * (ys[2] - ys[0]); - ys1[n-2] = 0.5 * (ys[n-1] - ys[n-3]); - ys1[n-1] = ys[n-1] - ys[n-2]; - for (int i = 2; i < n-2; i++) { - ys1[i]=((ys[i-2]-ys[i+2])+ 8.0*(ys[i+1]-ys[i-1]))/12.0; - } - for (int i = 0; i < n-1; i++) { - ys2[i]=3.0*(ys[i+1]-ys[i])-2.0*ys1[i]-ys1[i+1]; - ys3[i]=ys1[i]+ys1[i+1]-2.0*(ys[i+1]-ys[i]); - } - ys2[n-1]=0.0; - ys3[n-1]=0.0; - for (int i = 0; i < n; i++) { - ys4[i]=ys1[i]*rdx; - ys5[i]=2.0*ys2[i]*rdx; - ys6[i]=3.0*ys3[i]*rdx; - } - } - int size; - double xmin,xmax,xmaxsq,rdx,vmax; - double *ys, *ys1, *ys2, *ys3, *ys4, *ys5, *ys6; - double *xs; - }; -#endif diff --git a/src/tabular_function.h b/src/tabular_function.h index ec61c4bc4b..2931f296c8 100644 --- a/src/tabular_function.h +++ b/src/tabular_function.h @@ -18,11 +18,9 @@ namespace LAMMPS_NS { class TabularFunction { public: TabularFunction(); - TabularFunction(int, double, double); virtual ~TabularFunction(); void set_values(int, double, double, double *); - void set_values(int, double, double, double *, double); private: int size; @@ -34,6 +32,7 @@ class TabularFunction { void initialize(); public: + // used by pair style bop void value(double x, double &y, int ny, double &y1, int ny1) { double ps = (x - xmin) * rdx + 1.0; @@ -46,6 +45,7 @@ class TabularFunction { ys[ks - 1]; if (ny1) y1 = (ys6[ks - 1] * ps + ys5[ks - 1]) * ps + ys4[ks - 1]; } + // used by pair style polymorphic void value2(double x, double &y, int ny, double &y1, int ny1) { double ps = (x - xmin) * rdx;