include TabularFunction class
This commit is contained in:
Axel Kohlmeyer
279
src/tabular_function.cpp
Normal file
279
src/tabular_function.cpp
Normal file
@ -0,0 +1,279 @@
|
||||
/* ----------------------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
http://lammps.sandia.gov, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#include "tabular_function.h"
|
||||
|
||||
#include <cmath>
|
||||
#include <cstring>
|
||||
|
||||
using namespace LAMMPS_NS;
|
||||
|
||||
TabularFunction::TabularFunction()
|
||||
: size(0), xmin(0.0), xmax(0.0), xmaxsq(0.0), rdx(0.0), vmax(0.0),
|
||||
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;
|
||||
delete [] ys;
|
||||
delete [] ys1;
|
||||
delete [] ys2;
|
||||
delete [] ys3;
|
||||
delete [] ys4;
|
||||
delete [] ys5;
|
||||
delete [] ys6;
|
||||
}
|
||||
|
||||
void TabularFunction::set_values(int n, double x1, double x2, double *values)
|
||||
{
|
||||
reset_size(n);
|
||||
set_xrange(x1, x2);
|
||||
memcpy(ys, values, n*sizeof(double));
|
||||
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;
|
||||
xmax = x2;
|
||||
xmaxsq = x2*x2;
|
||||
}
|
||||
|
||||
void TabularFunction::reset_size(int n)
|
||||
{
|
||||
if (n != size) {
|
||||
size = n;
|
||||
delete [] xs;
|
||||
delete [] ys;
|
||||
delete [] ys1;
|
||||
delete [] ys2;
|
||||
delete [] ys3;
|
||||
delete [] ys4;
|
||||
delete [] ys5;
|
||||
delete [] ys6;
|
||||
xs = new double[n];
|
||||
ys = new double[n];
|
||||
ys1 = new double[n];
|
||||
ys2 = new double[n];
|
||||
ys3 = new double[n];
|
||||
ys4 = new double[n];
|
||||
ys5 = new double[n];
|
||||
ys6 = new double[n];
|
||||
}
|
||||
}
|
||||
|
||||
void TabularFunction::initialize()
|
||||
{
|
||||
int i;
|
||||
rdx = (xmax - xmin) / (size - 1.0);
|
||||
for (i = 0; i < size; 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[size - 2] = 0.5 * (ys[size - 1] - ys[size - 3]);
|
||||
ys1[size - 1] = ys[size - 1] - ys[size - 2];
|
||||
for (i = 2; i < size - 2; i++)
|
||||
ys1[i] = ((ys[i - 2] - ys[i + 2]) + 8.0 * (ys[i + 1] - ys[i - 1])) / 12.0;
|
||||
for (i = 0; i < size - 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[size - 1] = 0.0;
|
||||
ys3[size - 1] = 0.0;
|
||||
for (i = 0; i < size; i++) {
|
||||
ys4[i] = ys1[i] * rdx;
|
||||
ys5[i] = 2.0 * ys2[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
|
||||
68
src/tabular_function.h
Normal file
68
src/tabular_function.h
Normal file
@ -0,0 +1,68 @@
|
||||
/* -*- c++ -*- ----------------------------------------------------------
|
||||
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
|
||||
https://lammps.sandia.gov/, Sandia National Laboratories
|
||||
Steve Plimpton, sjplimp@sandia.gov
|
||||
|
||||
Copyright (2003) Sandia Corporation. Under the terms of Contract
|
||||
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
|
||||
certain rights in this software. This software is distributed under
|
||||
the GNU General Public License.
|
||||
|
||||
See the README file in the top-level LAMMPS directory.
|
||||
------------------------------------------------------------------------- */
|
||||
|
||||
#ifndef LMP_TABULAR_FUNCTION_H
|
||||
#define LMP_TABULAR_FUNCTION_H
|
||||
|
||||
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;
|
||||
double xmin, xmax, xmaxsq, rdx, vmax;
|
||||
double *xs, *ys, *ys1, *ys2, *ys3, *ys4, *ys5, *ys6;
|
||||
|
||||
void set_xrange(double x1, double x2);
|
||||
void reset_size(int);
|
||||
void initialize();
|
||||
|
||||
public:
|
||||
void value(double x, double &y, int ny, double &y1, int ny1)
|
||||
{
|
||||
double ps = (x - xmin) * rdx + 1.0;
|
||||
int ks = ps;
|
||||
if (ks > size - 1) ks = size - 1;
|
||||
ps = ps - ks;
|
||||
if (ps > 1.0) ps = 1.0;
|
||||
if (ny)
|
||||
y = ((ys3[ks - 1] * ps + ys2[ks - 1]) * ps + ys1[ks - 1]) * ps +
|
||||
ys[ks - 1];
|
||||
if (ny1) y1 = (ys6[ks - 1] * ps + ys5[ks - 1]) * ps + ys4[ks - 1];
|
||||
}
|
||||
void value2(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];
|
||||
}
|
||||
|
||||
double get_xmin() const { return xmin; }
|
||||
double get_xmax() const { return xmax; }
|
||||
double get_xmaxsq() const { return xmaxsq; }
|
||||
double get_rdx() const { return rdx; }
|
||||
double get_vmax() { return vmax; }
|
||||
};
|
||||
} // namespace LAMMPS_NS
|
||||
|
||||
#endif
|
||||
Reference in New Issue
Block a user