Geometry parameter calculation implemented
This commit is contained in:
Philipp Kloza
@ -13,6 +13,7 @@
|
||||
#include "integrate.h"
|
||||
#include "respa.h"
|
||||
#include "math_const.h"
|
||||
#include "math_extra.h"
|
||||
#include "memory.h"
|
||||
#include "error.h"
|
||||
#include "utils.h"
|
||||
@ -103,17 +104,22 @@ void PairMesoCNT::settings(int narg, char **arg)
|
||||
|
||||
void PairMesoCNT::coeff(int narg, char **arg)
|
||||
{
|
||||
if (narg != 9) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (narg != 10) error->all(FLERR,"Incorrect args for pair coefficients");
|
||||
if (!allocated) allocate();
|
||||
|
||||
sigma = force->numeric(FLERR,arg[2]);
|
||||
epsilon = force->numeric(FLERR,arg[3]);
|
||||
n_sigma = force->numeric(FLERR,arg[4]);
|
||||
n = force->inumeric(FLERR,arg[2]);
|
||||
sigma = force->numeric(FLERR,arg[3]);
|
||||
epsilon = force->numeric(FLERR,arg[4]);
|
||||
n_sigma = force->numeric(FLERR,arg[5]);
|
||||
|
||||
gamma_file = arg[5];
|
||||
uinf_file = arg[6];
|
||||
usemi_file = arg[7];
|
||||
phi_file = arg[8];
|
||||
gamma_file = arg[6];
|
||||
uinf_file = arg[7];
|
||||
usemi_file = arg[8];
|
||||
phi_file = arg[9];
|
||||
|
||||
radius = 1.421*3*n / MY_2PI;
|
||||
comega = 0.275*(1.0 - 1.0/(1.0 + 0.59*radius));
|
||||
ctheta = 0.35 + 0.0226*(radius - 6.785);
|
||||
|
||||
//Parse and bcast data
|
||||
int me;
|
||||
@ -150,14 +156,14 @@ double PairMesoCNT::spline(double x, double xstart, double dx,
|
||||
i = 0;
|
||||
// warn if argument below spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument below spline interval",cerror,ninput);
|
||||
sprintf(str,"Argument below spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
else if(i > coeff_size-1){
|
||||
i = coeff_size-1;
|
||||
// warn if argument above spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument above spline interval",cerror,ninput);
|
||||
sprintf(str,"Argument above spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
|
||||
@ -170,6 +176,66 @@ double PairMesoCNT::spline(double x, double xstart, double dx,
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairMesoCNT::spline(double x, double y, double *xstart, double ystart,
|
||||
double *dx, double dy, double ***coeff, int coeff_size)
|
||||
{
|
||||
int i = floor((y - ystart)/dy);
|
||||
if(i < 0){
|
||||
i = 0;
|
||||
// warn if argument below spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument below spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
else if(i > coeff_size-1){
|
||||
i = coeff_size-1;
|
||||
// warn if argument above spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument above spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
|
||||
double ylo = ystart + i*dy;
|
||||
double ybar = (y - ylo)/dy;
|
||||
|
||||
// compute coefficients in y
|
||||
|
||||
double a0, a1, a2, a3;
|
||||
double p0, p1, p2, p3;
|
||||
|
||||
p1 = spline(x,xstart[0],dx[0],coeff[0],coeff_size);
|
||||
p2 = spline(x,xstart[1],dx[1],coeff[1],coeff_size);
|
||||
|
||||
a0 = p1;
|
||||
|
||||
if(i == 0){
|
||||
p3 = spline(x,xstart[2],dx[2],coeff[2],coeff_size);
|
||||
|
||||
a1 = p2 - p1;
|
||||
a3 = 0.5*(p3 - 2*p2 + p1);
|
||||
a2 = -a3;
|
||||
}
|
||||
else if(i == coeff_size-2){
|
||||
p0 = spline(x,xstart[i-1],dx[i-1],coeff[i-1],coeff_size);
|
||||
|
||||
a1 = 0.5*(p2 - p0);
|
||||
a3 = 0.5*(p2 - 2*p1 + p0);
|
||||
a2 = -2*a3;
|
||||
}
|
||||
else{
|
||||
p0 = spline(x,xstart[i-1],dx[i-1],coeff[i-1],coeff_size);
|
||||
p3 = spline(x,xstart[2],dx[2],coeff[2],coeff_size);
|
||||
|
||||
a1 = 0.5*(p2 - p0);
|
||||
a2 = 0.5*(-p3 + 4*p2 - 5*p1 + 2*p0);
|
||||
a3 = 0.5*(p3 - 3*p2 + 3*p1 - p0);
|
||||
}
|
||||
|
||||
return a0 + ybar*(a1 + ybar*(a2 + a3*ybar));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairMesoCNT::dspline(double x, double xstart, double dx,
|
||||
double **coeff, int coeff_size)
|
||||
{
|
||||
@ -178,21 +244,141 @@ double PairMesoCNT::dspline(double x, double xstart, double dx,
|
||||
i = 0;
|
||||
// warn if argument below spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument below spline interval",cerror,ninput);
|
||||
sprintf(str,"Argument below spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
else if(i > coeff_size-1){
|
||||
i = coeff_size-1;
|
||||
// warn if argument above spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument above spline interval",cerror,ninput);
|
||||
sprintf(str,"Argument above spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
|
||||
double xlo = xstart + i*dx;
|
||||
double xbar = (x - xlo)/dx;
|
||||
|
||||
return coeff[i][1] + xbar*(2*coeff[i][2] + 3*xbar*coeff[i][3]);
|
||||
return (coeff[i][1] + xbar*(2*coeff[i][2] + 3*xbar*coeff[i][3])) / dx;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairMesoCNT::dxspline(double x, double y, double *xstart, double ystart,
|
||||
double *dx, double dy, double ***coeff, int coeff_size)
|
||||
{
|
||||
int i = floor((y - ystart)/dy);
|
||||
if(i < 0){
|
||||
i = 0;
|
||||
// warn if argument below spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument below spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
else if(i > coeff_size-1){
|
||||
i = coeff_size-1;
|
||||
// warn if argument above spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument above spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
|
||||
double ylo = ystart + i*dy;
|
||||
double ybar = (y - ylo)/dy;
|
||||
|
||||
// compute coefficients in y
|
||||
|
||||
double a0, a1, a2, a3;
|
||||
double p0, p1, p2, p3;
|
||||
|
||||
p1 = dspline(x,xstart[0],dx[0],coeff[0],coeff_size);
|
||||
p2 = dspline(x,xstart[1],dx[1],coeff[1],coeff_size);
|
||||
|
||||
a0 = p1;
|
||||
|
||||
if(i == 0){
|
||||
p3 = dspline(x,xstart[2],dx[2],coeff[2],coeff_size);
|
||||
|
||||
a1 = p2 - p1;
|
||||
a3 = 0.5*(p3 - 2*p2 + p1);
|
||||
a2 = -a3;
|
||||
}
|
||||
else if(i == coeff_size-2){
|
||||
p0 = dspline(x,xstart[i-1],dx[i-1],coeff[i-1],coeff_size);
|
||||
|
||||
a1 = 0.5*(p2 - p0);
|
||||
a3 = 0.5*(p2 - 2*p1 + p0);
|
||||
a2 = -2*a3;
|
||||
}
|
||||
else{
|
||||
p0 = dspline(x,xstart[i-1],dx[i-1],coeff[i-1],coeff_size);
|
||||
p3 = dspline(x,xstart[2],dx[2],coeff[2],coeff_size);
|
||||
|
||||
a1 = 0.5*(p2 - p0);
|
||||
a2 = 0.5*(-p3 + 4*p2 - 5*p1 + 2*p0);
|
||||
a3 = 0.5*(p3 - 3*p2 + 3*p1 - p0);
|
||||
}
|
||||
|
||||
return a0 + ybar*(a1 + ybar*(a2 + a3*ybar));
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairMesoCNT::dyspline(double x, double y, double *xstart, double ystart,
|
||||
double *dx, double dy, double ***coeff, int coeff_size)
|
||||
{
|
||||
int i = floor((y - ystart)/dy);
|
||||
if(i < 0){
|
||||
i = 0;
|
||||
// warn if argument below spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument below spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
else if(i > coeff_size-1){
|
||||
i = coeff_size-1;
|
||||
// warn if argument above spline range
|
||||
char str[128];
|
||||
sprintf(str,"Argument above spline interval");
|
||||
error->warning(FLERR,str);
|
||||
}
|
||||
|
||||
double ylo = ystart + i*dy;
|
||||
double ybar = (y - ylo)/dy;
|
||||
|
||||
// compute coefficients in y
|
||||
|
||||
double a0, a1, a2, a3;
|
||||
double p0, p1, p2, p3;
|
||||
|
||||
p1 = spline(x,xstart[0],dx[0],coeff[0],coeff_size);
|
||||
p2 = spline(x,xstart[1],dx[1],coeff[1],coeff_size);
|
||||
|
||||
a0 = p1;
|
||||
|
||||
if(i == 0){
|
||||
p3 = spline(x,xstart[2],dx[2],coeff[2],coeff_size);
|
||||
|
||||
a1 = p2 - p1;
|
||||
a3 = 0.5*(p3 - 2*p2 + p1);
|
||||
a2 = -a3;
|
||||
}
|
||||
else if(i == coeff_size-2){
|
||||
p0 = spline(x,xstart[i-1],dx[i-1],coeff[i-1],coeff_size);
|
||||
|
||||
a1 = 0.5*(p2 - p0);
|
||||
a3 = 0.5*(p2 - 2*p1 + p0);
|
||||
a2 = -2*a3;
|
||||
}
|
||||
else{
|
||||
p0 = spline(x,xstart[i-1],dx[i-1],coeff[i-1],coeff_size);
|
||||
p3 = spline(x,xstart[2],dx[2],coeff[2],coeff_size);
|
||||
|
||||
a1 = 0.5*(p2 - p0);
|
||||
a2 = 0.5*(-p3 + 4*p2 - 5*p1 + 2*p0);
|
||||
a3 = 0.5*(p3 - 3*p2 + 3*p1 - p0);
|
||||
}
|
||||
|
||||
return (a1 + ybar*(2*a2 + 3*a3*ybar)) / dy;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
@ -375,4 +561,150 @@ void PairMesoCNT::read_file(char *file, double **data,
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairMesoCNT::uinf(double h, double alpha, double xi1, double xi2)
|
||||
{
|
||||
double salpha = sin(alpha);
|
||||
double salphasq = salpha*salpha;
|
||||
if(salphasq < 1.0e-6){
|
||||
return (xi2 - xi1) * spline(h,start_uinf,del_uinf,uinf_coeff,pot_points);
|
||||
}
|
||||
else{
|
||||
double omega = 1.0 / (1.0 - comega*salphasq);
|
||||
double a = omega * salpha;
|
||||
double zeta1 = xi1 * a;
|
||||
double zeta2 = xi2 * a;
|
||||
|
||||
double phi1, phi2;
|
||||
if(zeta1 < 0) phi1 = -spline(-zeta1,h,startzeta_phi,starth_phi,
|
||||
delzeta_phi,delh_phi,phi_coeff,pot_points);
|
||||
else phi1 = spline(zeta1,h,startzeta_phi,starth_phi,
|
||||
delzeta_phi,delh_phi,phi_coeff,pot_points);
|
||||
|
||||
if(zeta2 < 0) phi2 = -spline(-zeta2,h,startzeta_phi,starth_phi,
|
||||
delzeta_phi,delh_phi,phi_coeff,pot_points);
|
||||
else phi2 = spline(zeta2,h,startzeta_phi,starth_phi,
|
||||
delzeta_phi,delh_phi,phi_coeff,pot_points);
|
||||
|
||||
double gamma_orth = spline(h,start_gamma,del_gamma,
|
||||
gamma_coeff,gamma_points);
|
||||
double gamma = 1.0 + (gamma_orth - 1.0)*salphasq;
|
||||
|
||||
return gamma * (phi2 - phi1) / a;
|
||||
}
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
double PairMesoCNT::usemi(double h, double alpha,
|
||||
double xi1, double xi2, double etaend)
|
||||
{
|
||||
double salpha = sin(alpha);
|
||||
double salphasq = salpha*salpha;
|
||||
double calpha = cos(alpha);
|
||||
double omega = 1.0 / (1.0 - comega*salphasq);
|
||||
double theta = 1.0 - ctheta*salphasq;
|
||||
|
||||
int points = 100;
|
||||
double delxi = (xi2 - xi1) / (points -1);
|
||||
|
||||
double g, hbar, etabar;
|
||||
double sum = 0;
|
||||
|
||||
// first and last term in sum
|
||||
|
||||
g = xi1 * omega * salpha;
|
||||
hbar = sqrt(h*h + g*g);
|
||||
etabar = xi1 * calpha - theta*etaend;
|
||||
sum += spline(hbar,etabar,starth_usemi,startxi_usemi,
|
||||
delh_usemi,delxi_usemi,usemi_coeff,pot_points);
|
||||
g = xi2 * omega * salpha;
|
||||
hbar = sqrt(h*h + g*g);
|
||||
etabar = xi2 * calpha - theta*etaend;
|
||||
sum += spline(hbar,etabar,starth_usemi,startxi_usemi,
|
||||
delh_usemi,delxi_usemi,usemi_coeff,pot_points);
|
||||
sum *= 0.5;
|
||||
|
||||
for(int i = 1; i < points-1; i++){
|
||||
double xibar = xi1 + i*delxi;
|
||||
g * xibar * omega * salpha;
|
||||
hbar = sqrt(h*h + g*g);
|
||||
etabar = xibar*calpha - theta*etaend;
|
||||
sum += spline(hbar,etabar,starth_usemi,startxi_usemi,
|
||||
delh_usemi,delxi_usemi,usemi_coeff,pot_points);
|
||||
}
|
||||
|
||||
double gamma_orth = spline(h,start_gamma,del_gamma,
|
||||
gamma_coeff,gamma_points);
|
||||
double gamma = 1.0 + (gamma_orth - 1.0)*salphasq;
|
||||
|
||||
return delxi * gamma * sum;
|
||||
}
|
||||
|
||||
/* ---------------------------------------------------------------------- */
|
||||
|
||||
void PairMesoCNT::geom(const double *r1, const double *r2,
|
||||
const double *p1, const double *p2, double *param)
|
||||
{
|
||||
using namespace MathExtra;
|
||||
double r[3], p[3], delr[3], l[3], m[3], rbar[3], pbar[3], delrbar[3];
|
||||
double psil[3], psim[3], dell_psim[3], delpsil_m[3];
|
||||
double delr1[3], delr2[3], delp1[3], delp2[3];
|
||||
double ex[3], ey[3];
|
||||
double psi, frac, taur, taup;
|
||||
double h, alpha, xi1, xi2, eta1, eta2;
|
||||
|
||||
add3(r1,r2,r);
|
||||
scale3(0.5,r);
|
||||
add3(p1,p2,p);
|
||||
scale3(0.5,p);
|
||||
|
||||
sub3(p,r,delr);
|
||||
|
||||
sub3(r2,r1,l);
|
||||
normalize3(l,l);
|
||||
sub3(p2,p1,m);
|
||||
normalize3(m,m);
|
||||
|
||||
psi = dot3(l,m);
|
||||
frac = 1.0 / (1.0 - psi*psi);
|
||||
|
||||
copy3(l,psil);
|
||||
scale3(psi,psil);
|
||||
copy3(m,psim);
|
||||
scale3(psi,psim);
|
||||
|
||||
sub3(l,psim,dell_psim);
|
||||
sub3(psil,m,delpsil_m);
|
||||
taur = dot3(delr,dell_psim) * frac;
|
||||
taup = dot3(delr,delpsil_m) * frac;
|
||||
|
||||
scaleadd3(taur,l,r,rbar);
|
||||
scaleadd3(taup,m,p,pbar);
|
||||
sub3(pbar,rbar,delrbar);
|
||||
|
||||
h = len3(delrbar);
|
||||
|
||||
copy3(delrbar,ex);
|
||||
scale3(1/h,ex);
|
||||
cross3(l,ex,ey);
|
||||
|
||||
if(dot3(m,ey) < 0) alpha = acos(psi);
|
||||
else alpha = MY_2PI - acos(psi);
|
||||
|
||||
sub3(r1,rbar,delr1);
|
||||
sub3(r2,rbar,delr2);
|
||||
xi1 = dot3(delr1,l);
|
||||
xi2 = dot3(delr2,l);
|
||||
|
||||
sub3(p1,pbar,delp1);
|
||||
sub3(p2,pbar,delp2);
|
||||
eta1 = dot3(delp1,m);
|
||||
eta2 = dot3(delp2,m);
|
||||
|
||||
param[0] = h;
|
||||
param[1] = alpha;
|
||||
param[2] = xi1;
|
||||
param[3] = xi2;
|
||||
param[4] = eta1;
|
||||
param[5] = eta2;
|
||||
}
|
||||
|
||||
@ -18,23 +18,15 @@ class PairMesoCNT : public Pair {
|
||||
void compute(int, int);
|
||||
void settings(int, char **);
|
||||
void coeff(int, char **);
|
||||
void init_style();
|
||||
double init_one(int, int);
|
||||
void write_restart(FILE *);
|
||||
void read_restart(FILE *);
|
||||
|
||||
void write_restart_settings(FILE *);
|
||||
void read_restart_settings(FILE *);
|
||||
void write_data(FILE *);
|
||||
void write_data_all(FILE *);
|
||||
double single(int, int, int, int, double, double, double, double &);
|
||||
void *extract(const char *, int &);
|
||||
|
||||
//void init_style();
|
||||
//double init_one(int, int);
|
||||
|
||||
protected:
|
||||
int n, gamma_points, pot_points;
|
||||
double sigma, epsilon, n_sigma, radius;
|
||||
double sigma, epsilon, n_sigma, radius, comega, ctheta;
|
||||
double start_gamma, start_uinf, startxi_usemi, starth_phi;
|
||||
double del_gamma, del_uinf, delxi_usemi, delh_phi;
|
||||
double *starth_usemi, *startzeta_phi;
|
||||
double *delh_usemi, *delzeta_phi;
|
||||
double *gamma_data, *uinf_data;
|
||||
double **usemi_data, **phi_data;
|
||||
@ -45,16 +37,25 @@ class PairMesoCNT : public Pair {
|
||||
void allocate();
|
||||
|
||||
double spline(double, double, double, double **, int);
|
||||
double spline(double, double, double *, double, double *, double ***, int);
|
||||
double spline(double, double, double *, double, double *,
|
||||
double, double ***, int);
|
||||
double dspline(double, double, double, double **, int);
|
||||
double dxspline(double, double, double *, double, double *, double ***, int);
|
||||
double dyspline(double, double, double *, double, double *, double ***, int);
|
||||
double dxspline(double, double, double *, double, double *,
|
||||
double, double ***, int);
|
||||
double dyspline(double, double, double *, double, double *,
|
||||
double, double ***, int);
|
||||
|
||||
void spline_coeff(double *, double **, int);
|
||||
void spline_coeff(double **, double ***, int);
|
||||
|
||||
void read_file(char *, double *, double *, int);
|
||||
void read_file(char *, double **, double *, double *, int);
|
||||
|
||||
double uinf(double, double, double, double);
|
||||
double usemi(double, double, double, double, double);
|
||||
|
||||
void geom(const double *, const double *, const double *, const double *,
|
||||
double *);
|
||||
};
|
||||
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user