// // CDDL HEADER START // // The contents of this file are subject to the terms of the Common Development // and Distribution License Version 1.0 (the "License"). // // You can obtain a copy of the license at // http://www.opensource.org/licenses/CDDL-1.0. See the License for the // specific language governing permissions and limitations under the License. // // When distributing Covered Code, include this CDDL HEADER in each file and // include the License file in a prominent location with the name LICENSE.CDDL. // If applicable, add the following below this CDDL HEADER, with the fields // enclosed by brackets "[]" replaced with your own identifying information: // // Portions Copyright (c) [yyyy] [name of copyright owner]. All rights reserved. // // CDDL HEADER END // // // Copyright (c) 2014, Regents of the University of Minnesota. // All rights reserved. // // Contributors: // Ryan S. Elliott // Mingjian Wen // #include "EAM_QuinticClampedSpline.hpp" #include "EAM_Implementation.hpp" //****************************************************************************** void EAM_Implementation::SplineInterpolate(double const * const dat, double const delta, int const n, double * const coe) { // setup convenient pointers (spline) into the coefficients (coe) array double ** const spline = new double *[n]; // deleted at end of function for (int i = 0; i < n; ++i) { spline[i] = &coe[i * NUMBER_SPLINE_COEFF]; } // B.C.s end first and second derivatives four point finite difference spline[0][F_LINEAR] = (-11 * dat[0] + 18 * dat[1] - 9 * dat[2] + 2 * dat[3]) / 6; spline[0][F_QUADRATIC] = (2 * dat[0] - 5 * dat[1] + 4 * dat[2] - dat[3]) / 2; spline[n - 1][F_LINEAR] = (-2 * dat[n - 4] + 9 * dat[n - 3] - 18 * dat[n - 2] + 11 * dat[n - 1]) / 6; spline[n - 1][F_QUADRATIC] = (-dat[n - 4] + 4 * dat[n - 3] - 5 * dat[n - 2] + 2 * dat[n - 1]) / 2; // local variables double * a = new double[2 * (n - 2)]; double * b = new double[2 * (n - 2)]; double * c = new double[2 * (n - 2)]; double * d = new double[2 * (n - 2)]; double * e = new double[2 * (n - 2)]; double * f = new double[2 * (n - 2)]; double * g = new double[2 * (n - 2)]; double * h = new double[2 * (n - 2)]; double * xn = new double[2 * (n - 2)]; double xmult; // put the matrix values into the variables for (int i = 0; i < n - 2; ++i) { d[2 * i] = 6; d[2 * i + 1] = 16; c[2 * i] = 0; c[2 * i + 1] = -4; e[2 * i] = 0; e[2 * i + 1] = -2; b[2 * i] = -1; b[2 * i + 1] = 7; f[2 * i] = -1; f[2 * i + 1] = 7; a[2 * i] = 2; a[2 * i + 1] = 0; g[2 * i] = 4; g[2 * i + 1] = 0; h[2 * i] = 10 * (dat[i + 2] - 2 * dat[i + 1] + dat[i]); h[2 * i + 1] = 15 * (dat[i + 2] - dat[i]); } h[0] = h[0] + spline[0][F_QUADRATIC] + 4 * spline[0][F_LINEAR]; h[1] = h[1] - 2 * spline[0][F_QUADRATIC] - 7 * spline[0][F_LINEAR]; h[2 * (n - 2) - 2] = h[2 * (n - 2) - 2] + spline[n - 1][F_QUADRATIC] - 4 * spline[n - 1][F_LINEAR]; h[2 * (n - 2) - 1] = h[2 * (n - 2) - 1] + 2 * spline[n - 1][F_QUADRATIC] - 7 * spline[n - 1][F_LINEAR]; // gaussian forward elimination for (int i = 1; i < 2 * (n - 2) - 2; ++i) { xmult = c[i - 1] / d[i - 1]; d[i] = d[i] - xmult * e[i - 1]; e[i] = e[i] - xmult * f[i - 1]; f[i] = f[i] - xmult * g[i - 1]; h[i] = h[i] - xmult * h[i - 1]; xmult = b[i - 1] / d[i - 1]; c[i] = c[i] - xmult * e[i - 1]; d[i + 1] = d[i + 1] - xmult * f[i - 1]; e[i + 1] = e[i + 1] - xmult * g[i - 1]; h[i + 1] = h[i + 1] - xmult * h[i - 1]; xmult = a[i - 1] / d[i - 1]; b[i] = b[i] - xmult * e[i - 1]; c[i + 1] = c[i + 1] - xmult * f[i - 1]; d[i + 2] = d[i + 2] - xmult * g[i - 1]; h[i + 2] = h[i + 2] - xmult * h[i - 1]; } // elimination of the last two lines xmult = c[2 * (n - 2) - 3] / d[2 * (n - 2) - 3]; d[2 * (n - 2) - 2] = d[2 * (n - 2) - 2] - xmult * e[2 * (n - 2) - 3]; e[2 * (n - 2) - 2] = e[2 * (n - 2) - 2] - xmult * f[2 * (n - 2) - 3]; h[2 * (n - 2) - 2] = h[2 * (n - 2) - 2] - xmult * h[2 * (n - 2) - 3]; xmult = b[2 * (n - 2) - 3] / d[2 * (n - 2) - 3]; c[2 * (n - 2) - 2] = c[2 * (n - 2) - 2] - xmult * e[2 * (n - 2) - 3]; d[2 * (n - 2) - 1] = d[2 * (n - 2) - 1] - xmult * f[2 * (n - 2) - 3]; h[2 * (n - 2) - 1] = h[2 * (n - 2) - 1] - xmult * h[2 * (n - 2) - 3]; // eliminaiton of the last line xmult = c[2 * (n - 2) - 2] / d[2 * (n - 2) - 2]; d[2 * (n - 2) - 1] = d[2 * (n - 2) - 1] - xmult * e[2 * (n - 2) - 2]; h[2 * (n - 2) - 1] = h[2 * (n - 2) - 1] - xmult * h[2 * (n - 2) - 2]; // gaussian backward substitution xn[2 * (n - 2) - 1] = h[2 * (n - 2) - 1] / d[2 * (n - 2) - 1]; xn[2 * (n - 2) - 2] = (h[2 * (n - 2) - 2] - e[2 * (n - 2) - 2] * xn[2 * (n - 2) - 1]) / d[2 * (n - 2) - 2]; xn[2 * (n - 2) - 3] = (h[2 * (n - 2) - 3] - e[2 * (n - 2) - 3] * xn[2 * (n - 2) - 2] - f[2 * (n - 2) - 3] * xn[2 * (n - 2) - 1]) / d[2 * (n - 2) - 3]; for (int i = 2 * (n - 2) - 4; i >= 0; --i) { xn[i] = (h[i] - e[i] * xn[i + 1] - f[i] * xn[i + 2] - g[i] * xn[i + 3]) / d[i]; } // put the coefficinets into coe for (int i = 1; i < n - 1; i++) { spline[i][F_QUADRATIC] = xn[2 * i - 2]; spline[i][F_LINEAR] = xn[2 * i - 1]; } for (int m = 0; m < n - 1; ++m) { spline[m][F_CONSTANT] = dat[m]; spline[m][F_CUBIC] = 10 * (dat[m + 1] - dat[m]) - 6 * spline[m][F_LINEAR] - 4 * spline[m + 1][F_LINEAR] - 3 * spline[m][F_QUADRATIC] + spline[m + 1][F_QUADRATIC]; spline[m][F_QUARTIC] = -15 * (dat[m + 1] - dat[m]) + 8 * spline[m][F_LINEAR] + 7 * spline[m + 1][F_LINEAR] + 3 * spline[m][F_QUADRATIC] - 2 * spline[m + 1][F_QUADRATIC]; spline[m][F_QUINTIC] = 6 * (dat[m + 1] - dat[m]) - 3 * spline[m][F_LINEAR] - 3 * spline[m + 1][F_LINEAR] - spline[m][F_QUADRATIC] + spline[m + 1][F_QUADRATIC]; } for (int m = 0; m < n - 1; ++m) { spline[m][DF_CONSTANT] = spline[m][F_LINEAR] / delta; spline[m][DF_LINEAR] = 2.0 * spline[m][F_QUADRATIC] / delta; spline[m][DF_QUADRATIC] = 3.0 * spline[m][F_CUBIC] / delta; spline[m][DF_CUBIC] = 4.0 * spline[m][F_QUARTIC] / delta; spline[m][DF_QUARTIC] = 5.0 * spline[m][F_QUINTIC] / delta; } for (int m = 0; m < n - 1; ++m) { spline[m][D2F_CONSTANT] = spline[m][DF_LINEAR] / delta; spline[m][D2F_LINEAR] = 2.0 * spline[m][DF_QUADRATIC] / delta; spline[m][D2F_QUADRATIC] = 3.0 * spline[m][DF_CUBIC] / delta; spline[m][D2F_CUBIC] = 4.0 * spline[m][DF_QUARTIC] / delta; } delete[] a; delete[] b; delete[] c; delete[] d; delete[] e; delete[] f; delete[] g; delete[] h; delete[] xn; delete[] spline; }