ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
e58aeebcce3e15fc315c34354ddb4eae5f6dc5cf
lammps/src/USER-SMD/fix_smd_wall_surface.cpp

489 lines
16 KiB
C++
Raw Blame History

/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
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.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Mike Parks (SNL), Ezwanur Rahman, J.T. Foster (UTSA)
------------------------------------------------------------------------- */
#include "fix_smd_wall_surface.h"
#include <mpi.h>
#include <cstring>
#include <cstdlib>
#include <Eigen/Eigen>
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "error.h"
#include "atom_vec.h"
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace Eigen;
using namespace std;
#define DELTA 16384
#define EPSILON 1.0e-6
/* ---------------------------------------------------------------------- */
FixSMDWallSurface::FixSMDWallSurface(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
restart_global = 0;
restart_peratom = 0;
first = 1;
//atom->add_callback(0);
//atom->add_callback(1);
if (narg != 6)
error->all(FLERR, "Illegal number of arguments for fix smd/wall_surface");
filename = strdup(arg[3]);
wall_particle_type = force->inumeric(FLERR, arg[4]);
wall_molecule_id = force->inumeric(FLERR, arg[5]);
if (wall_molecule_id < 65535) {
error->one(FLERR, "wall molcule id must be >= 65535\n");
}
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix smd/wall_surface reads trianglulated surface from file: %s\n", filename);
printf("fix smd/wall_surface has particle type %d \n", wall_particle_type);
printf("fix smd/wall_surface has molecule id %d \n", wall_molecule_id);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
}
/* ---------------------------------------------------------------------- */
FixSMDWallSurface::~FixSMDWallSurface() {
free(filename);
filename = NULL;
// unregister this fix so atom class doesn't invoke it any more
//atom->delete_callback(id, 0);
//atom->delete_callback(id, 1);
}
/* ---------------------------------------------------------------------- */
int FixSMDWallSurface::setmask() {
int mask = 0;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDWallSurface::init() {
if (!first)
return;
}
/* ----------------------------------------------------------------------
For minimization: setup as with dynamics
------------------------------------------------------------------------- */
void FixSMDWallSurface::min_setup(int vflag) {
setup(vflag);
}
/* ----------------------------------------------------------------------
create initial list of neighbor partners via call to neighbor->build()
must be done in setup (not init) since fix init comes before neigh init
------------------------------------------------------------------------- */
void FixSMDWallSurface::setup(int /*vflag*/) {
if (!first)
return;
first = 0;
// set bounds for my proc
// if periodic and I am lo/hi proc, adjust bounds by EPSILON
// insures all data atoms will be owned even with round-off
int triclinic = domain->triclinic;
double epsilon[3];
if (triclinic)
epsilon[0] = epsilon[1] = epsilon[2] = EPSILON;
else {
epsilon[0] = domain->prd[0] * EPSILON;
epsilon[1] = domain->prd[1] * EPSILON;
epsilon[2] = domain->prd[2] * EPSILON;
}
if (triclinic == 0) {
sublo[0] = domain->sublo[0];
subhi[0] = domain->subhi[0];
sublo[1] = domain->sublo[1];
subhi[1] = domain->subhi[1];
sublo[2] = domain->sublo[2];
subhi[2] = domain->subhi[2];
} else {
sublo[0] = domain->sublo_lamda[0];
subhi[0] = domain->subhi_lamda[0];
sublo[1] = domain->sublo_lamda[1];
subhi[1] = domain->subhi_lamda[1];
sublo[2] = domain->sublo_lamda[2];
subhi[2] = domain->subhi_lamda[2];
}
if (comm->layout != Comm::LAYOUT_TILED) {
if (domain->xperiodic) {
if (comm->myloc[0] == 0)
sublo[0] -= epsilon[0];
if (comm->myloc[0] == comm->procgrid[0] - 1)
subhi[0] += epsilon[0];
}
if (domain->yperiodic) {
if (comm->myloc[1] == 0)
sublo[1] -= epsilon[1];
if (comm->myloc[1] == comm->procgrid[1] - 1)
subhi[1] += epsilon[1];
}
if (domain->zperiodic) {
if (comm->myloc[2] == 0)
sublo[2] -= epsilon[2];
if (comm->myloc[2] == comm->procgrid[2] - 1)
subhi[2] += epsilon[2];
}
} else {
if (domain->xperiodic) {
if (comm->mysplit[0][0] == 0.0)
sublo[0] -= epsilon[0];
if (comm->mysplit[0][1] == 1.0)
subhi[0] += epsilon[0];
}
if (domain->yperiodic) {
if (comm->mysplit[1][0] == 0.0)
sublo[1] -= epsilon[1];
if (comm->mysplit[1][1] == 1.0)
subhi[1] += epsilon[1];
}
if (domain->zperiodic) {
if (comm->mysplit[2][0] == 0.0)
sublo[2] -= epsilon[2];
if (comm->mysplit[2][1] == 1.0)
subhi[2] += epsilon[2];
}
}
read_triangles(0);
}
/* ----------------------------------------------------------------------
function to determine number of values in a text line
------------------------------------------------------------------------- */
int FixSMDWallSurface::count_words(const char *line) {
int n = strlen(line) + 1;
char *copy;
memory->create(copy, n, "atom:copy");
strcpy(copy, line);
char *ptr;
if ((ptr = strchr(copy, '#')))
*ptr = '\0';
if (strtok(copy, " \t\n\r\f") == NULL) {
memory->destroy(copy);
return 0;
}
n = 1;
while (strtok(NULL, " \t\n\r\f"))
n++;
memory->destroy(copy);
return n;
}
/* ----------------------------------------------------------------------
size of atom nlocal's restart data
------------------------------------------------------------------------- */
void FixSMDWallSurface::read_triangles(int pass) {
double coord[3];
int nlocal_previous = atom->nlocal;
int ilocal = nlocal_previous;
int m;
int me;
bigint natoms_previous = atom->natoms;
Vector3d *vert;
vert = new Vector3d[3];
Vector3d normal, center;
FILE *fp = fopen(filename, "r");
if (fp == NULL) {
char str[128];
snprintf(str,128, "Cannot open file %s", filename);
error->one(FLERR, str);
}
MPI_Comm_rank(world, &me);
if (me == 0) {
if (screen) {
if (pass == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
fprintf(screen, " scanning triangle pairs ...\n");
} else {
fprintf(screen, " reading triangle pairs ...\n");
}
}
if (logfile) {
if (pass == 0) {
fprintf(logfile, " scanning triangle pairs ...\n");
} else {
fprintf(logfile, " reading triangle pairs ...\n");
}
}
}
char line[256];
char *retpointer;
char **values;
int nwords;
// read STL solid name
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
error->one(FLERR,"error reading number of triangle pairs");
}
nwords = count_words(line);
if (nwords < 1) {
error->one(FLERR,"first line of file is incorrect");
}
// values = new char*[nwords];
// values[0] = strtok(line, " \t\n\r\f");
// if (values[0] == NULL)
// error->all(FLERR, "Incorrect atom format in data file");
// for (m = 1; m < nwords; m++) {
// values[m] = strtok(NULL, " \t\n\r\f");
// if (values[m] == NULL)
// error->all(FLERR, "Incorrect atom format in data file");
// }
// delete[] values;
//
// if (comm->me == 0) {
// cout << "STL file contains solid body with name: " << values[1] << endl;
// }
// iterate over STL facets util end of body is reached
while (fgets(line, sizeof(line), fp)) { // read a line, should be the facet line
// evaluate facet line
nwords = count_words(line);
if (nwords != 5) {
//sprintf(str, "found end solid line");
//error->message(FLERR, str);
break;
} else {
// should be facet line
}
values = new char*[nwords];
values[0] = strtok(line, " \t\n\r\f");
if (values[0] == NULL)
error->all(FLERR, "Incorrect atom format in data file");
for (m = 1; m < nwords; m++) {
values[m] = strtok(NULL, " \t\n\r\f");
if (values[m] == NULL)
error->all(FLERR, "Incorrect atom format in data file");
}
normal << force->numeric(FLERR, values[2]), force->numeric(FLERR, values[3]), force->numeric(FLERR, values[4]);
//cout << "normal is " << normal << endl;
delete[] values;
// read outer loop line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
error->one(FLERR, "error reading outer loop");
}
nwords = count_words(line);
if (nwords != 2) {
error->one(FLERR,"error reading outer loop");
}
// read vertex lines
for (int k = 0; k < 3; k++) {
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
error->one(FLERR,"error reading vertex line");
}
nwords = count_words(line);
if (nwords != 4) {
error->one(FLERR,"error reading vertex line");
}
values = new char*[nwords];
values[0] = strtok(line, " \t\n\r\f");
if (values[0] == NULL)
error->all(FLERR,"Incorrect vertex line");
for (m = 1; m < nwords; m++) {
values[m] = strtok(NULL, " \t\n\r\f");
if (values[m] == NULL)
error->all(FLERR, "Incorrect vertex line");
}
vert[k] << force->numeric(FLERR, values[1]), force->numeric(FLERR, values[2]), force->numeric(FLERR, values[3]);
//cout << "vertex is " << vert[k] << endl;
//printf("%s %s %s\n", values[1], values[2], values[3]);
delete[] values;
//exit(1);
}
// read end loop line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
error->one(FLERR, "error reading endloop");
}
nwords = count_words(line);
if (nwords != 1) {
error->one(FLERR,"error reading endloop");
}
// read end facet line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
error->one(FLERR,"error reading endfacet");
}
nwords = count_words(line);
if (nwords != 1) {
error->one(FLERR,"error reading endfacet");
}
// now we have a normal and three vertices ... proceed with adding triangle
center = (vert[0] + vert[1] + vert[2]) / 3.0;
// cout << "center is " << center << endl;
double r1 = (center - vert[0]).norm();
double r2 = (center - vert[1]).norm();
double r3 = (center - vert[2]).norm();
double r = MAX(r1, r2);
r = MAX(r, r3);
/*
* if atom/molecule is in my subbox, create it
* ... use x0 to hold triangle normal.
* ... use smd_data_9 to hold the three vertices
* ... use x to hold triangle center
* ... radius is the mmaximal distance from triangle center to all vertices
*/
// printf("coord: %f %f %f\n", coord[0], coord[1], coord[2]);
// printf("sublo: %f %f %f\n", sublo[0], sublo[1], sublo[2]);
// printf("subhi: %f %f %f\n", subhi[0], subhi[1], subhi[2]);
//printf("ilocal = %d\n", ilocal);
if (center(0) >= sublo[0] && center(0) < subhi[0] && center(1) >= sublo[1] && center(1) < subhi[1] && center(2) >= sublo[2]
&& center(2) < subhi[2]) {
//printf("******* KERATIN nlocal=%d ***\n", nlocal);
coord[0] = center(0);
coord[1] = center(1);
coord[2] = center(2);
atom->avec->create_atom(wall_particle_type, coord);
/*
* need to initialize pointers to atom vec arrays here, because they could have changed
* due to calling grow() in create_atoms() above;
*/
tagint *mol = atom->molecule;
int *type = atom->type;
double *radius = atom->radius;
double *contact_radius = atom->contact_radius;
double **smd_data_9 = atom->smd_data_9;
double **x0 = atom->x0;
radius[ilocal] = r; //ilocal;
contact_radius[ilocal] = r; //ilocal;
mol[ilocal] = wall_molecule_id;
type[ilocal] = wall_particle_type;
x0[ilocal][0] = normal(0);
x0[ilocal][1] = normal(1);
x0[ilocal][2] = normal(2);
smd_data_9[ilocal][0] = vert[0](0);
smd_data_9[ilocal][1] = vert[0](1);
smd_data_9[ilocal][2] = vert[0](2);
smd_data_9[ilocal][3] = vert[1](0);
smd_data_9[ilocal][4] = vert[1](1);
smd_data_9[ilocal][5] = vert[1](2);
smd_data_9[ilocal][6] = vert[2](0);
smd_data_9[ilocal][7] = vert[2](1);
smd_data_9[ilocal][8] = vert[2](2);
ilocal++;
}
}
// set new total # of atoms and error check
bigint nblocal = atom->nlocal;
MPI_Allreduce(&nblocal, &atom->natoms, 1, MPI_LMP_BIGINT, MPI_SUM, world);
if (atom->natoms < 0 || atom->natoms >= MAXBIGINT)
error->all(FLERR, "Too many total atoms");
// add IDs for newly created atoms
// check that atom IDs are valid
if (atom->tag_enable)
atom->tag_extend();
atom->tag_check();
// create global mapping of atoms
// zero nghost in case are adding new atoms to existing atoms
if (atom->map_style) {
atom->nghost = 0;
atom->map_init();
atom->map_set();
}
// print status
if (comm->me == 0) {
if (screen) {
printf("... fix smd/wall_surface finished reading triangulated surface\n");
fprintf(screen, "fix smd/wall_surface created " BIGINT_FORMAT " atoms\n", atom->natoms - natoms_previous);
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
if (logfile) {
fprintf(logfile, "... fix smd/wall_surface finished reading triangulated surface\n");
fprintf(logfile, "fix smd/wall_surface created " BIGINT_FORMAT " atoms\n", atom->natoms - natoms_previous);
fprintf(logfile, ">>========>>========>>========>>========>>========>>========>>========>>========\n");
}
}
delete[] vert;
fclose(fp);
}
Reference in New Issue View Git Blame Copy Permalink