ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
3ea061279597bb352aed0946746ebaba72cff01b
lammps/src/MANIFOLD/manifold_thylakoid.cpp
Axel Kohlmeyer 2132b1d904 update developer reference text
2022-10-24 11:08:26 -04:00

496 lines
13 KiB
C++
Raw Blame History

// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
LAMMPS development team: developers@lammps.org
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 "manifold_thylakoid.h"
#include "manifold_thylakoid_shared.h"
#include "comm.h"
#include "domain.h" // For some checks regarding the simulation box.
#include "error.h"
#include <cmath>
using namespace LAMMPS_NS;
using namespace user_manifold;
manifold_thylakoid::manifold_thylakoid( LAMMPS *lmp, int /*narg*/, char ** /*arg*/)
: manifold(lmp)
{
// You can NOT depend on proper construction of the domains in
// the constructor, because the params are set by the function that
// calls the factory-construction method. Instead, the constructing
// fix should call post_param_init();
}
manifold_thylakoid::~manifold_thylakoid()
{
for (auto & part : parts) {
delete part;
}
}
void manifold_thylakoid::post_param_init()
{
// Set coefficients:
lT = 3; // Radius of cylinder edges of stacks
LT = 15; // Size of faces of cylinder stacks
pad = 3; // Padding to prevent improper look-ups
wB = 3.0;
LB = 10.0;
lB = 3.0;
wB = params[0];
LB = params[1];
lB = params[2];
init_domains();
}
double manifold_thylakoid::g( const double *x )
{
int err = 0;
std::size_t idx;
thyla_part *p = get_thyla_part(x,&err,&idx);
if (err) error->one(FLERR,"Error getting thyla_part for x = ({}, {}, {})",x[0],x[1],x[2]);
double con_val = p->g(x);
if (std::isfinite(con_val)) {
return con_val;
} else {
error->one(FLERR,"Error, thyla_part of type {} returned {} as constraint val!", p->type, con_val);
return 0;
}
}
void manifold_thylakoid::n( const double *x, double *n )
{
int err = 0;
std::size_t idx;
thyla_part *p = get_thyla_part(x,&err,&idx);
if (err)
error->one(FLERR,"Error getting thyla_part for x = ({}, {}, {})",x[0],x[1],x[2]);
p->n(x,n);
if (std::isfinite(n[0]) && std::isfinite(n[1]) && std::isfinite(n[2])) {
return;
} else
error->one(FLERR,"thyla_part of type {} returned ({},{},{}) as gradient!",
p->type, n[0], n[1], n[2]);
}
thyla_part *manifold_thylakoid::get_thyla_part( const double *x, int * /*err_flag*/, std::size_t *idx )
{
for (std::size_t i = 0; i < parts.size(); ++i) {
thyla_part *p = parts[i];
if (is_in_domain(p,x)) {
if (idx != nullptr) *idx = i;
return p;
}
}
error->one(FLERR,"Could not find thyla_part for x = ({},{},{})", x[0],x[1],x[2]);
return nullptr;
}
void manifold_thylakoid::init_domains()
{
if (wB + 2*lB > LT)
error->one(FLERR,"LT = {} not large enough to accommodate bridge with "
"wB = {} and lB = {}! {} > {}\n", LT, wB, lB, wB + 2*lB, LT);
// Determine some constant coordinates:
x0 = -( 0.5*LB + lB + lT + LT + lT + pad);
y0 = -( 0.5*LT + lT + pad );
z0 = -15;
if (x0 < domain->boxlo[0])
error->one(FLERR,"Thylakoid expects xlo of at most {:.8f}, but found {:.8f}", x0, domain->boxlo[0]);
if (y0 < domain->boxlo[1])
error->one(FLERR,"Thylakoid expects ylo of at most {:.8f}, but found {:.8f}",y0, domain->boxlo[1]);
// Add some padding to prevent improper lookups.
z0 -= pad;
x1 = -x0;
y1 = -y0;
z1 = -z0;
Lx = x1 - x0;
Ly = y1 - y0;
Lz = z1 - z0;
if (x1 > domain->boxhi[0])
error->one(FLERR,"Expected xhi larger than current box has: {:.8f} > {:.8f}", x1, domain->boxhi[0]);
if (y1 > domain->boxhi[1])
error->one(FLERR,"Expected yhi larger than current box has: {:.8f} > {:.8f}", y1, domain->boxhi[1]);
// Create and add the manifold parts to the array.
thyla_part *p;
// Determine coordinates of domain boundaries and centres of "mass":
thyla_part_geom cllb, cllt, clrb, clrt; // Left thylakoid cylinder parts
thyla_part_geom pll, plb, plt, plr; // Left thylakoid plane parts
thyla_part_geom crlb, crlt, crrb, crrt; // Right thylakoid cylinder parts
thyla_part_geom prl, prb, prt, prr; // Right thylakoid plane parts
thyla_part_geom bl, br, bc; // Bridge left, right connectors and cylinder.
// The bridge is three parts.
// 1. A connector between the right face of the left grana and a cylinder
// 2. A connector between the left face of the right grana and a cylinder
// 3. The aforementioned cylinder.
// 1.:
// Args: pt, X0, R0, R, x0, y0, z0, sign
double rB = 0.5*wB;
double Reff = rB + lB;
double safety_fac = 1.0;
bl.pt[0] = -0.5*LB;
bl.pt[1] = 0;
bl.pt[2] = 0;
bl.lo[0] = bl.pt[0] - safety_fac*lB;
bl.lo[1] = -(1.0 + safety_fac) * Reff;
bl.lo[2] = -(1.0 + safety_fac) * Reff;
bl.hi[0] = bl.pt[0];
bl.hi[1] = (1.0 + safety_fac) * Reff;
bl.hi[2] = (1.0 + safety_fac) * Reff;
// double X0, double R0, double R, double s,
p = make_cyl_to_plane_part(bl.pt[0], rB, lB, -1, bl.pt);
set_domain(p, bl.lo, bl.hi);
parts.push_back(p);
// 2.:
br.pt[0] = 0.5*LB;
br.pt[1] = 0;
br.pt[2] = 0;
br.lo[0] = br.pt[0];
br.lo[1] = -(1.0 + safety_fac) * Reff;
br.lo[2] = -(1.0 + safety_fac) * Reff;
br.hi[0] = br.pt[0] + safety_fac*lB;
br.hi[1] = (1.0 + safety_fac) * Reff;
br.hi[2] = (1.0 + safety_fac) * Reff;
// double X0, double R0, double R, double s,
p = make_cyl_to_plane_part(br.pt[0], rB, lB, 1, br.pt);
set_domain(p, br.lo, br.hi);
parts.push_back(p);
// 3.:
// Cylinder in between:
bc.pt[0] = 0;
bc.pt[1] = 0;
bc.pt[2] = 0;
bc.lo[0] = bl.pt[0];
bc.lo[1] = -Reff;
bc.lo[2] = -Reff;
bc.hi[0] = br.pt[0];
bc.hi[1] = Reff;
bc.hi[2] = Reff;
p = make_cyl_part( 0, 1, 1, bc.pt, rB );
set_domain( p, bc.lo, bc.hi );
parts.push_back(p);
// The stack on the left:
cllb.lo[0] = x0;
cllb.lo[1] = y0;
cllb.lo[2] = z0;
cllb.pt[0] = x0 + pad + lT;
cllb.pt[1] = y0 + pad + lT;
cllb.pt[2] = 0;
cllb.hi[0] = cllb.pt[0];
cllb.hi[1] = cllb.pt[1];
cllb.hi[2] = z1;
p = make_cyl_part(1,1,0,cllb.pt,lT);
set_domain(p, cllb.lo, cllb.hi);
parts.push_back(p);
// left left top cylinder
cllt = cllb;
cllt.lo[1] = y1 - pad - lT;
cllt.hi[1] = y1;
cllt.pt[1] = cllb.pt[1] + LT;
p = make_cyl_part(1,1,0,cllt.pt,lT);
set_domain(p, cllt.lo, cllt.hi);
parts.push_back(p);
// left right bottom cylinder
clrb = cllb;
clrb.pt[0] += LT;
clrb.lo[0] = clrb.pt[0];
clrb.hi[0] = clrb.lo[0] + lT + lB;
p = make_cyl_part(1,1,0,clrb.pt,lT);
set_domain(p, clrb.lo, clrb.hi);
parts.push_back(p);
// left right top cylinder
clrt = clrb;
clrt.pt[1] += LT;
clrt.lo[1] = y1 - pad - lT;
clrt.hi[1] = y1;
p = make_cyl_part(1,1,0,clrt.pt,lT);
set_domain(p, clrt.lo, clrt.hi);
parts.push_back(p);
// left left plane
pll.pt[0] = x0 + pad;
pll.pt[1] = 0;
pll.pt[2] = 0;
pll.lo[0] = x0;
pll.lo[1] = cllb.pt[1];
pll.lo[2] = z0;
pll.hi[0] = pll.lo[0] + pad + lT;
pll.hi[1] = pll.lo[1] + LT;
pll.hi[2] = z1;
p = make_plane_part(1,0,0,pll.pt);
set_domain(p, pll.lo, pll.hi);
parts.push_back(p);
// left bottom plane
plb.pt[0] = x0 + pad + lT + 0.5*LT;
plb.pt[1] = y0 + pad;
plb.pt[2] = 0;
plb.lo[0] = x0 + pad + lT;
plb.lo[1] = y0;
plb.lo[2] = z0;
plb.hi[0] = plb.lo[0] + LT;
plb.hi[1] = plb.lo[1] + pad + lT;
plb.hi[2] = z1;
p = make_plane_part(0,1,0,plb.pt);
set_domain(p, plb.lo, plb.hi);
parts.push_back(p);
// left top plane
plt = plb;
plt.lo[1] = cllb.pt[1] + LT;
plt.hi[1] = y1;
plt.pt[1] = y1 - pad;
p = make_plane_part(0,1,0,plt.pt);
set_domain(p, plt.lo, plt.hi);
parts.push_back(p);
// left right plane
plr = pll;
plr.lo[0] = bl.lo[0] - 0.5;
plr.lo[1] = y0 - pad;
plr.hi[0] = bl.lo[0] + 0.5;
plr.hi[1] = y1 + pad;
plr.pt[0] = bl.pt[0] - lB;
plr.pt[1] = 0.0;
plr.pt[2] = 0.0;
plr.hi[2] = z1 + pad;
plr.lo[2] = z0 - pad;
p = make_plane_part(1,0,0,plr.pt);
set_domain(p, plr.lo, plr.hi);
parts.push_back(p);
// Check if this plane lines up with bl:
if (fabs(plr.pt[0] - bl.pt[0] + lB) > 1e-8)
error->one(FLERR,"Origins of plane left right and bridge left misaligned! {} != {}!\n",
plr.pt[0], bl.pt[0] - lB );
// Now, for the right stack, you can mirror the other...
// To mirror them you need to invert lo[0] and hi[0] and flip their sign.
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crlb, &cllb);
p = make_cyl_part(1,1,0,crlb.pt,lT);
set_domain(p, crlb.lo, crlb.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crlt, &cllt);
p = make_cyl_part(1,1,0,crlt.pt,lT);
set_domain(p, crlt.lo, crlt.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crrb, &clrb);
p = make_cyl_part(1,1,0,crrb.pt,lT);
set_domain(p, crrb.lo, crrb.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &crrt, &clrt);
p = make_cyl_part(1,1,0,crrt.pt,lT);
set_domain(p, crrt.lo, crrt.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prl , &pll );
p = make_plane_part(1,0,0,prl.pt);
set_domain(p, prl.lo, prl.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prb , &plb );
p = make_plane_part(0,1,0,prb.pt);
set_domain(p, prb.lo, prb.hi);
parts.push_back(p);
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prt , &plt );
p = make_plane_part(0,1,0,prt.pt);
set_domain(p, prt.lo, prt.hi);
parts.push_back(p);
// Careful, this one is wrongly named.
thyla_part_geom::mirror(thyla_part_geom::DIR_X, &prr, &plr);
p = make_plane_part(1,0,0,prr.pt);
set_domain(p, prr.lo, prr.hi);
parts.push_back(p);
if (fabs(prr.pt[0] - br.pt[0] - lB) > 1e-8)
error->one(FLERR,"Origins of plane left right and bridge left misaligned! {} != {}!\n",
prr.pt[0], br.pt[0] + lB);
}
void manifold_thylakoid::set_domain( thyla_part *p, const std::vector<double> &lo,
const std::vector<double> &hi )
{
if (lo[0] >= hi[0]) {
error->one(FLERR,"xlo >= xhi ({:.8f} >= {:.8f})",lo[0],hi[0]);
} else if (lo[1] >= hi[1]) {
error->one(FLERR,"ylo >= yhi ({:.8f} >= {:.8f})",lo[1],hi[1]);
} else if (lo[2] >= hi[2]) {
error->one(FLERR,"zlo >= zhi ({:.8f} >= {:.8f})",lo[2],hi[2]);
}
p->xlo = lo[0];
p->ylo = lo[1];
p->zlo = lo[2];
p->xhi = hi[0];
p->yhi = hi[1];
p->zhi = hi[2];
}
int manifold_thylakoid::is_in_domain( thyla_part *part, const double *x )
{
bool domain_ok = (x[0] >= part->xlo) && (x[0] <= part->xhi) &&
(x[1] >= part->ylo) && (x[1] <= part->yhi) &&
(x[2] >= part->zlo) && (x[2] <= part->zhi);
if (!domain_ok) return false;
// From here on out, domain is ok.
if (part->type == thyla_part::THYLA_TYPE_CYL_TO_PLANE) {
double R0 = part->params[1];
double R = part->params[2];
double y = x[1];
double z = x[2];
double dist2 = y*y + z*z;
double RR = R0+R;
double RR2 = RR*RR;
if (dist2 < RR2) {
return true;
} else {
// Domain was ok, but radius not.
return false;
}
} else {
return true;
}
}
thyla_part *manifold_thylakoid::make_plane_part (double a, double b, double c,
const std::vector<double> &pt )
{
double args[7];
args[0] = a;
args[1] = b;
args[2] = c;
args[3] = pt[0];
args[4] = pt[1];
args[5] = pt[2];
auto p = new thyla_part(thyla_part::THYLA_TYPE_PLANE,args,0,0,0,0,0,0);
return p;
}
thyla_part *manifold_thylakoid::make_cyl_part (double a, double b, double c,
const std::vector<double> &pt, double R)
{
double args[7];
args[0] = a;
args[1] = b;
args[2] = c;
args[3] = pt[0];
args[4] = pt[1];
args[5] = pt[2];
args[6] = R;
auto p = new thyla_part(thyla_part::THYLA_TYPE_CYL,args,0,0,0,0,0,0);
return p;
}
thyla_part *manifold_thylakoid::make_sphere_part(const std::vector<double> &pt, double R)
{
double args[7];
args[0] = R;
args[1] = pt[0];
args[2] = pt[1];
args[3] = pt[2];
auto p = new thyla_part(thyla_part::THYLA_TYPE_SPHERE,args,0,0,0,0,0,0);
return p;
}
thyla_part *manifold_thylakoid::make_cyl_to_plane_part(double X0, double R0, double R,
double s, const std::vector<double> &pt )
{
double args[7];
args[0] = X0;
args[1] = R0;
args[2] = R;
args[3] = pt[0];
args[4] = pt[1];
args[5] = pt[2];
args[6] = s;
auto p = new thyla_part(thyla_part::THYLA_TYPE_CYL_TO_PLANE,args,0,0,0,0,0,0);
return p;
}
Reference in New Issue View Git Blame Copy Permalink