Homec4science
Diffusion lammps (efficient_neuronet) Browse

lammps/src/USER-MANIFOLD/manifold_thylakoid_shared.cpp
49804d97059fefficient_neuronet

Last Change
Raw File

manifold_thylakoid_shared.cpp
View Options

#include "manifold_thylakoid_shared.h"
#include <math.h>
using namespace LAMMPS_NS;
using namespace user_manifold;
thyla_part::thyla_part( int type, double *args, double xlo, double ylo, double zlo,
double xhi, double yhi, double zhi )
: type(type), xlo(xlo), xhi(xhi),
ylo(ylo), yhi(yhi), zlo(zlo), zhi(zhi)
{
switch(type){
case THYLA_TYPE_PLANE: // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0
params[0] = args[0]; // a
params[1] = args[1]; // b
params[2] = args[2]; // c
params[3] = args[3]; // x0
params[4] = args[4]; // y0
params[5] = args[5]; // z0
break;
case THYLA_TYPE_SPHERE: // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0
params[0] = args[0]; // R
params[1] = args[1]; // x0
params[2] = args[2]; // y0
params[3] = args[3]; // z0
break;
case THYLA_TYPE_CYL: // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0
params[0] = args[0]; // a
params[1] = args[1]; // b
params[2] = args[2]; // c
params[3] = args[3]; // x0
params[4] = args[4]; // y0
params[5] = args[5]; // z0
params[6] = args[6]; // R
if( (args[0] != 0.0) && (args[1] != 0.0) && (args[2] != 0.0) ){
err_flag = -1;
return;
}
// The others should be 1.
if( (args[0] != 1.0) && (args[0] != 0.0) &&
(args[1] != 1.0) && (args[1] != 0.0) &&
(args[2] != 1.0) && (args[2] != 0.0) ){
err_flag = -1;
}
break;
case THYLA_TYPE_CYL_TO_PLANE:
/*
* Funky bit that connects a cylinder to a plane.
* It is what you get by rotating the equation
* r(x) = R0 + R*( 1 - sqrt( 1 - ( ( X0 - x ) /R )^2 ) ) around the x-axis.
* I kid you not.
*
* The shape is symmetric so you have to set whether it is the "left" or
* "right" end by truncating it properly. It is designed to run from
* X0 to X0 + R if "right" or X0 - R to X0 if "left".
*
* As params it takes X0, R0, R, and a sign that determines whether it is
* "left" (args[3] < 0) or "right" (args[3] > 0).
*
* The args are: X0, R0, R, x0, y0, z0, sign.
*/
params[0] = args[0];
params[1] = args[1];
params[2] = args[2];
params[3] = args[3];
params[4] = args[4];
params[5] = args[5];
params[6] = args[6];
break;
default:
err_flag = -1;
}
x0 = (type == THYLA_TYPE_SPHERE) ? params[1] : params[3];
y0 = (type == THYLA_TYPE_SPHERE) ? params[2] : params[4];
z0 = (type == THYLA_TYPE_SPHERE) ? params[3] : params[5];
}
thyla_part::~thyla_part()
{}
double thyla_part::g(const double *x)
{
switch(type){
case THYLA_TYPE_PLANE:{ // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0
double a = params[0];
double b = params[1];
double c = params[2];
double dx = x[0] - params[3];
double dy = x[1] - params[4];
double dz = x[2] - params[5];
return a*dx + b*dy + c*dz;
break;
}
case THYLA_TYPE_SPHERE:{ // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0
double R2 = params[0]*params[0];
double dx = x[0] - params[1];
double dy = x[1] - params[2];
double dz = x[2] - params[3];
return dx*dx + dy*dy + dz*dz - R2;
break;
}
case THYLA_TYPE_CYL:{ // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0
double a = params[0];
double b = params[1];
double c = params[2];
double X0 = params[3];
double Y0 = params[4];
double Z0 = params[5];
double R = params[6];
double dx = x[0] - X0;
double dy = x[1] - Y0;
double dz = x[2] - Z0;
return a*dx*dx + b*dy*dy + c*dz*dz - R*R;
break;
}
case THYLA_TYPE_CYL_TO_PLANE:{
double X0 = params[0];
double R0 = params[1];
double R = params[2];
// Determine the centre of the sphere.
double dx = (x[0] - X0);
double dyz = sqrt( x[1]*x[1] + x[2]*x[2] );
double rdyz = dyz - (R0 + R);
double norm = 1.0 / (2.0*R);
// Maybe sign is important here...
double g = norm*(dx*dx + rdyz*rdyz - R*R);
return g;
}
default:
err_flag = -1;
return 0; // Mostly to get rid of compiler werrors.
break;
}
}
void thyla_part::n( const double *x, double *n )
{
switch(type){
case THYLA_TYPE_PLANE:{ // a*(x-x0) + b*(y-y0) + c*(z-z0) = 0
double a = params[0];
double b = params[1];
double c = params[2];
n[0] = a;
n[1] = b;
n[2] = c;
break;
}
case THYLA_TYPE_SPHERE:{ // (x-x0)^2 + (y-y0)^2 + (z-z0)^2 - R^2 = 0
double dx = x[0] - params[1];
double dy = x[1] - params[2];
double dz = x[2] - params[3];
n[0] = 2*dx;
n[1] = 2*dy;
n[2] = 2*dz;
break;
}
case THYLA_TYPE_CYL:{ // a*(x-x0)^2 + b*(y-y0)^2 + c*(z-z0)^2 - R^2 = 0
double a = params[0];
double b = params[1];
double c = params[2];
double X0 = params[3];
double Y0 = params[4];
double Z0 = params[5];
double dx = x[0] - X0;
double dy = x[1] - Y0;
double dz = x[2] - Z0;
n[0] = 2*a*dx;
n[1] = 2*b*dy;
n[2] = 2*c*dz;
break;
}
case THYLA_TYPE_CYL_TO_PLANE:{
double X0 = params[0];
double R0 = params[1];
double R = params[2];
double s = (params[6] > 0.0) ? 1.0 : -1.0;
// Determine the centre of the sphere.
double dx = s*(x[0] - X0);
double ryz = sqrt( x[1]*x[1] + x[2]*x[2] );
// Maybe sign is important here...
// Normalize g and n so that the normal is continuous:
double norm = 1.0 / (2.0 * R);
n[0] = s*2*dx*norm;
double const_part = 1.0 - (R0 + R) / ryz;
n[1] = 2*x[1]*const_part*norm;
n[2] = 2*x[2]*const_part*norm;
break;
}
default:
err_flag = -1;
break;
}
}
void thyla_part_geom::mirror( unsigned int axis, thyla_part_geom *m,
const thyla_part_geom *o )
{
// Since dir is already the index of the array this is really simple:
m->lo[0] = o->lo[0];
m->lo[1] = o->lo[1];
m->lo[2] = o->lo[2];
m->pt[0] = o->pt[0];
m->pt[1] = o->pt[1];
m->pt[2] = o->pt[2];
m->hi[0] = o->hi[0];
m->hi[1] = o->hi[1];
m->hi[2] = o->hi[2];
m->lo[axis] = -o->hi[axis];
m->hi[axis] = -o->lo[axis];
m->pt[axis] = -o->pt[axis];
}
c4science · Help