From 4f9f87c8ab4dc0efa3cad28be77dc1cd3bdb7b1d Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Tue, 4 May 2021 00:24:12 -0400 Subject: [PATCH] include TabularFunction class --- src/tabular_function.cpp | 279 +++++++++++++++++++++++++++++++++++++++ src/tabular_function.h | 68 ++++++++++ 2 files changed, 347 insertions(+) create mode 100644 src/tabular_function.cpp create mode 100644 src/tabular_function.h diff --git a/src/tabular_function.cpp b/src/tabular_function.cpp new file mode 100644 index 0000000000..3c4922b92b --- /dev/null +++ b/src/tabular_function.cpp @@ -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 +#include + +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 diff --git a/src/tabular_function.h b/src/tabular_function.h new file mode 100644 index 0000000000..ec61c4bc4b --- /dev/null +++ b/src/tabular_function.h @@ -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