ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

simplify remove unused code

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-05-04 22:26:02 -04:00
parent 09684b776f
commit 8285e71fd9
2 changed files with 6 additions and 182 deletions
Download Patch File Download Diff File Expand all files Collapse all files

184
src/tabular_function.cpp
View File

@ -23,13 +23,6 @@ TabularFunction::TabularFunction()
xs(nullptr), ys(nullptr), ys1(nullptr), ys2(nullptr), ys3(nullptr), xs(nullptr), ys(nullptr), ys1(nullptr), ys2(nullptr), ys3(nullptr),
ys4(nullptr), ys5(nullptr), ys6(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() TabularFunction::~TabularFunction()
{ {
delete [] xs; delete [] xs;
@ -50,42 +43,6 @@ void TabularFunction::set_values(int n, double x1, double x2, double *values)
initialize(); 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) void TabularFunction::set_xrange(double x1, double x2)
{ {
xmin = x1; xmin = x1;
@ -120,6 +77,10 @@ void TabularFunction::initialize()
{ {
int i; int i;
rdx = (xmax - xmin) / (size - 1.0); 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; for (i = 0; i < size; i++) xs[i] = xmin + i * rdx;
rdx = 1.0 / rdx; rdx = 1.0 / rdx;
ys1[0] = ys[1] - ys[0]; ys1[0] = ys[1] - ys[0];
@ -140,140 +101,3 @@ void TabularFunction::initialize()
ys6[i] = 3.0 * ys3[i] * rdx; 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

4
src/tabular_function.h
View File

@ -18,11 +18,9 @@ namespace LAMMPS_NS {
class TabularFunction { class TabularFunction {
public: public:
TabularFunction(); TabularFunction();
TabularFunction(int, double, double);
virtual ~TabularFunction(); virtual ~TabularFunction();
void set_values(int, double, double, double *); void set_values(int, double, double, double *);
void set_values(int, double, double, double *, double);
private: private:
int size; int size;
@ -34,6 +32,7 @@ class TabularFunction {
void initialize(); void initialize();
public: public:
// used by pair style bop
void value(double x, double &y, int ny, double &y1, int ny1) void value(double x, double &y, int ny, double &y1, int ny1)
{ {
double ps = (x - xmin) * rdx + 1.0; double ps = (x - xmin) * rdx + 1.0;
@ -46,6 +45,7 @@ class TabularFunction {
ys[ks - 1]; ys[ks - 1];
if (ny1) y1 = (ys6[ks - 1] * ps + ys5[ks - 1]) * ps + ys4[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) void value2(double x, double &y, int ny, double &y1, int ny1)
{ {
double ps = (x - xmin) * rdx; double ps = (x - xmin) * rdx;