1 /* ----------------------------------------------------------------------
2 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
3 https://www.lammps.org/ Sandia National Laboratories
4 Steve Plimpton, sjplimp@sandia.gov
5
6 Copyright (2003) Sandia Corporation. Under the terms of Contract
7 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
8 certain rights in this software. This software is distributed under
9 the GNU General Public License.
10
11 See the README file in the top-level LAMMPS directory.
12 ------------------------------------------------------------------------- */
13
14 #include "tabular_function.h"
15
16 #include <cmath>
17 #include <cstring>
18
19 using namespace LAMMPS_NS;
20
TabularFunction()21 TabularFunction::TabularFunction() :
22 size(0), xmin(0.0), xmax(0.0), xmaxsq(0.0), rdx(0.0), vmax(0.0), xs(nullptr), ys(nullptr),
23 ys1(nullptr), ys2(nullptr), ys3(nullptr), ys4(nullptr), ys5(nullptr), ys6(nullptr)
24 {
25 }
26
~TabularFunction()27 TabularFunction::~TabularFunction()
28 {
29 delete[] xs;
30 delete[] ys;
31 delete[] ys1;
32 delete[] ys2;
33 delete[] ys3;
34 delete[] ys4;
35 delete[] ys5;
36 delete[] ys6;
37 }
38
set_values(int n,double x1,double x2,double * values)39 void TabularFunction::set_values(int n, double x1, double x2, double *values)
40 {
41 reset_size(n);
42 set_xrange(x1, x2);
43 memcpy(ys, values, n * sizeof(double));
44 initialize();
45 }
46
set_xrange(double x1,double x2)47 void TabularFunction::set_xrange(double x1, double x2)
48 {
49 xmin = x1;
50 xmax = x2;
51 xmaxsq = x2 * x2;
52 }
53
reset_size(int n)54 void TabularFunction::reset_size(int n)
55 {
56 if (n != size) {
57 size = n;
58 delete[] xs;
59 delete[] ys;
60 delete[] ys1;
61 delete[] ys2;
62 delete[] ys3;
63 delete[] ys4;
64 delete[] ys5;
65 delete[] ys6;
66 xs = new double[n];
67 ys = new double[n];
68 ys1 = new double[n];
69 ys2 = new double[n];
70 ys3 = new double[n];
71 ys4 = new double[n];
72 ys5 = new double[n];
73 ys6 = new double[n];
74 }
75 }
76
initialize()77 void TabularFunction::initialize()
78 {
79 int i;
80 rdx = (xmax - xmin) / (size - 1.0);
81 vmax = 0.0;
82 for (i = 0; i < size; i++)
83 if (fabs(ys[i]) > vmax) vmax = fabs(ys[i]);
84 for (i = 0; i < size; i++) xs[i] = xmin + i * rdx;
85 rdx = 1.0 / rdx;
86 ys1[0] = ys[1] - ys[0];
87 ys1[1] = 0.5 * (ys[2] - ys[0]);
88 ys1[size - 2] = 0.5 * (ys[size - 1] - ys[size - 3]);
89 ys1[size - 1] = ys[size - 1] - ys[size - 2];
90 for (i = 2; i < size - 2; i++)
91 ys1[i] = ((ys[i - 2] - ys[i + 2]) + 8.0 * (ys[i + 1] - ys[i - 1])) / 12.0;
92 for (i = 0; i < size - 1; i++) {
93 ys2[i] = 3.0 * (ys[i + 1] - ys[i]) - 2.0 * ys1[i] - ys1[i + 1];
94 ys3[i] = ys1[i] + ys1[i + 1] - 2.0 * (ys[i + 1] - ys[i]);
95 }
96 ys2[size - 1] = 0.0;
97 ys3[size - 1] = 0.0;
98 for (i = 0; i < size; i++) {
99 ys4[i] = ys1[i] * rdx;
100 ys5[i] = 2.0 * ys2[i] * rdx;
101 ys6[i] = 3.0 * ys3[i] * rdx;
102 }
103 }
104