ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge pull request #3717 from evoyiatzis/master

Browse Source 
Block regions with bounds defined by equal-style variables
This commit is contained in:
Axel Kohlmeyer
2023-04-03 21:31:17 -04:00
committed by GitHub
parent c9b6fad6a0 751dbd08b6
commit cbf7f66a6f
3 changed files with 211 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/src/region.rst
View File

@ -18,6 +18,7 @@ Syntax
*delete* = no args
*block* args = xlo xhi ylo yhi zlo zhi
xlo,xhi,ylo,yhi,zlo,zhi = bounds of block in all dimensions (distance units)
xlo,xhi,ylo,yhi,zlo,zhi can be a variable
*cone* args = dim c1 c2 radlo radhi lo hi
dim = *x* or *y* or *z* = axis of cone
c1,c2 = coords of cone axis in other 2 dimensions (distance units)

204
src/region_block.cpp
View File

@ -16,68 +16,114 @@
#include "domain.h"
#include "error.h"
#include "input.h"
#include "math_extra.h"
#include "variable.h"
#include <cstring>
using namespace LAMMPS_NS;
enum{CONSTANT,VARIABLE};
#define BIG 1.0e20
/* ---------------------------------------------------------------------- */
RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) :
Region(lmp, narg, arg), xlostr(nullptr), xhistr(nullptr), ylostr(nullptr), yhistr(nullptr), zlostr(nullptr), zhistr(nullptr)
{
options(narg-8,&arg[8]);
xlostyle = CONSTANT;
if (strcmp(arg[2],"INF") == 0 || strcmp(arg[2],"EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[2],"INF") == 0) xlo = -BIG;
else if (domain->triclinic == 0) xlo = domain->boxlo[0];
else xlo = domain->boxlo_bound[0];
} else if (utils::strmatch(arg[2],"^v_")) {
xlostr = utils::strdup(arg[2]+2);
xlo = 0.0;
xlostyle = VARIABLE;
varshape = 1;
} else xlo = xscale*utils::numeric(FLERR,arg[2],false,lmp);
xhistyle = CONSTANT;
if (strcmp(arg[3],"INF") == 0 || strcmp(arg[3],"EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[3],"INF") == 0) xhi = BIG;
else if (domain->triclinic == 0) xhi = domain->boxhi[0];
else xhi = domain->boxhi_bound[0];
} else if (utils::strmatch(arg[3],"^v_")) {
xhistr = utils::strdup(arg[3]+2);
xhi = 0.0;
xhistyle = VARIABLE;
varshape = 1;
} else xhi = xscale*utils::numeric(FLERR,arg[3],false,lmp);
ylostyle = CONSTANT;
if (strcmp(arg[4],"INF") == 0 || strcmp(arg[4],"EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[4],"INF") == 0) ylo = -BIG;
else if (domain->triclinic == 0) ylo = domain->boxlo[1];
else ylo = domain->boxlo_bound[1];
} else if (utils::strmatch(arg[4],"^v_")) {
ylostr = utils::strdup(arg[4]+2);
ylo = 0.0;
ylostyle = VARIABLE;
varshape = 1;
} else ylo = yscale*utils::numeric(FLERR,arg[4],false,lmp);
yhistyle = CONSTANT;
if (strcmp(arg[5],"INF") == 0 || strcmp(arg[5],"EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[5],"INF") == 0) yhi = BIG;
else if (domain->triclinic == 0) yhi = domain->boxhi[1];
else yhi = domain->boxhi_bound[1];
} else if (utils::strmatch(arg[5],"^v_")) {
yhistr = utils::strdup(arg[5]+2);
yhi = 0.0;
yhistyle = VARIABLE;
varshape = 1;
} else yhi = yscale*utils::numeric(FLERR,arg[5],false,lmp);
zlostyle = CONSTANT;
if (strcmp(arg[6],"INF") == 0 || strcmp(arg[6],"EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[6],"INF") == 0) zlo = -BIG;
else if (domain->triclinic == 0) zlo = domain->boxlo[2];
else zlo = domain->boxlo_bound[2];
} else if (utils::strmatch(arg[6],"^v_")) {
zlostr = utils::strdup(arg[6]+2);
zlo = 0.0;
zlostyle = VARIABLE;
varshape = 1;
} else zlo = zscale*utils::numeric(FLERR,arg[6],false,lmp);
zhistyle = CONSTANT;
if (strcmp(arg[7],"INF") == 0 || strcmp(arg[7],"EDGE") == 0) {
if (domain->box_exist == 0)
error->all(FLERR,"Cannot use region INF or EDGE when box does not exist");
if (strcmp(arg[7],"INF") == 0) zhi = BIG;
else if (domain->triclinic == 0) zhi = domain->boxhi[2];
else zhi = domain->boxhi_bound[2];
} else if (utils::strmatch(arg[7],"^v_")) {
zhistr = utils::strdup(arg[7]+2);
zhi = 0.0;
zhistyle = VARIABLE;
varshape = 1;
} else zhi = zscale*utils::numeric(FLERR,arg[7],false,lmp);
if (varshape) {
variable_check();
RegBlock::shape_update();
}
// error check
if (xlo > xhi) error->all(FLERR,"Illegal region block xlo: {} >= xhi: {}", xlo, xhi);
@ -189,10 +235,23 @@ RegBlock::RegBlock(LAMMPS *lmp, int narg, char **arg) : Region(lmp, narg, arg)
RegBlock::~RegBlock()
{
if (copymode) return;
delete [] xlostr;
delete [] xhistr;
delete [] ylostr;
delete [] yhistr;
delete [] zlostr;
delete [] zhistr;
delete [] contact;
}
/* ---------------------------------------------------------------------- */
void RegBlock::init()
{
Region::init();
if (varshape) variable_check();
}
/* ----------------------------------------------------------------------
inside = 1 if x,y,z is inside or on surface
inside = 0 if x,y,z is outside and not on surface
@ -340,6 +399,147 @@ int RegBlock::surface_exterior(double *x, double cutoff)
return 0;
}
/* ----------------------------------------------------------------------
change region shape via variable evaluation
------------------------------------------------------------------------- */
void RegBlock::shape_update() // addition
{
if (xlostyle == VARIABLE)
xlo = xscale * input->variable->compute_equal(xlovar);
if (xhistyle == VARIABLE)
xhi = xscale * input->variable->compute_equal(xhivar);
if (ylostyle == VARIABLE)
ylo = yscale * input->variable->compute_equal(ylovar);
if (yhistyle == VARIABLE)
yhi = yscale * input->variable->compute_equal(yhivar);
if (zlostyle == VARIABLE)
zlo = zscale * input->variable->compute_equal(zlovar);
if (zhistyle == VARIABLE)
zhi = zscale * input->variable->compute_equal(zhivar);
if (xlo > xhi || ylo > yhi || zlo > zhi)
error->one(FLERR,"Variable evaluation in region gave bad value");
// face[0]
corners[0][0][0] = xlo;
corners[0][0][1] = ylo;
corners[0][0][2] = zlo;
corners[0][1][0] = xlo;
corners[0][1][1] = ylo;
corners[0][1][2] = zhi;
corners[0][2][0] = xlo;
corners[0][2][1] = yhi;
corners[0][2][2] = zhi;
corners[0][3][0] = xlo;
corners[0][3][1] = yhi;
corners[0][3][2] = zlo;
// face[1]
corners[1][0][0] = xhi;
corners[1][0][1] = ylo;
corners[1][0][2] = zlo;
corners[1][1][0] = xhi;
corners[1][1][1] = ylo;
corners[1][1][2] = zhi;
corners[1][2][0] = xhi;
corners[1][2][1] = yhi;
corners[1][2][2] = zhi;
corners[1][3][0] = xhi;
corners[1][3][1] = yhi;
corners[1][3][2] = zlo;
// face[2]
MathExtra::copy3(corners[0][0], corners[2][0]);
MathExtra::copy3(corners[1][0], corners[2][1]);
MathExtra::copy3(corners[1][1], corners[2][2]);
MathExtra::copy3(corners[0][1], corners[2][3]);
// face[3]
MathExtra::copy3(corners[0][3], corners[3][0]);
MathExtra::copy3(corners[0][2], corners[3][1]);
MathExtra::copy3(corners[1][2], corners[3][2]);
MathExtra::copy3(corners[1][3], corners[3][3]);
// face[4]
MathExtra::copy3(corners[0][0], corners[4][0]);
MathExtra::copy3(corners[0][3], corners[4][1]);
MathExtra::copy3(corners[1][3], corners[4][2]);
MathExtra::copy3(corners[1][0], corners[4][3]);
// face[5]
MathExtra::copy3(corners[0][1], corners[5][0]);
MathExtra::copy3(corners[1][1], corners[5][1]);
MathExtra::copy3(corners[1][2], corners[5][2]);
MathExtra::copy3(corners[0][2], corners[5][3]);
}
/* ----------------------------------------------------------------------
error check on existence of variable
------------------------------------------------------------------------- */
void RegBlock::variable_check() // addition
{
if (xlostyle == VARIABLE) {
xlovar = input->variable->find(xlostr);
if (xlovar < 0)
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(xlovar))
error->all(FLERR,"Variable for region block is invalid style");
}
if (xhistyle == VARIABLE) {
xhivar = input->variable->find(xhistr);
if (xhivar < 0)
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(xhivar))
error->all(FLERR,"Variable for region block is invalid style");
}
if (ylostyle == VARIABLE) {
ylovar = input->variable->find(ylostr);
if (ylovar < 0)
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(ylovar))
error->all(FLERR,"Variable for region block is invalid style");
}
if (yhistyle == VARIABLE) {
yhivar = input->variable->find(yhistr);
if (yhivar < 0)
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(yhivar))
error->all(FLERR,"Variable for region block is invalid style");
}
if (zlostyle == VARIABLE) {
zlovar = input->variable->find(zlostr);
if (zlovar < 0)
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(zlovar))
error->all(FLERR,"Variable for region block is invalid style");
}
if (zhistyle == VARIABLE) {
zhivar = input->variable->find(zhistr);
if (zhivar < 0)
error->all(FLERR,"Variable name for region block does not exist");
if (!input->variable->equalstyle(zhivar))
error->all(FLERR,"Variable for region block is invalid style");
}
}
/*------------------------------------------------------------------------
return distance to closest point on surface I of block region
store closest point in xc,yc,zc

8
src/region_block.h
View File

@ -30,17 +30,25 @@ class RegBlock : public Region {
public:
RegBlock(class LAMMPS *, int, char **);
~RegBlock() override;
void init() override;
int inside(double, double, double) override;
int surface_interior(double *, double) override;
int surface_exterior(double *, double) override;
void shape_update() override;
protected:
double xlo, xhi, ylo, yhi, zlo, zhi;
double corners[6][4][3];
double face[6][3];
int xlostyle, xlovar, xhistyle, xhivar;
int ylostyle, ylovar, yhistyle, yhivar;
int zlostyle, zlovar, zhistyle, zhivar;
char *xlostr, *ylostr, *zlostr;
char *xhistr, *yhistr, *zhistr;
double find_closest_point(int, double *, double &, double &, double &);
int inside_face(double *, int);
void variable_check();
};
} // namespace LAMMPS_NS