Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F85159916
fix_smd_wall_surface.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Fri, Sep 27, 04:48
Size
13 KB
Mime Type
text/x-c
Expires
Sun, Sep 29, 04:48 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
21113630
Attached To
rLAMMPS lammps
fix_smd_wall_surface.cpp
View Options
/* ----------------------------------------------------------------------
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 <math.h>
#include "fix_smd_wall_surface.h"
#include "atom.h"
#include "domain.h"
#include "force.h"
#include "comm.h"
#include "update.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "pair.h"
#include "lattice.h"
#include "memory.h"
#include "error.h"
#include <Eigen/Eigen>
#include <stdio.h>
#include "atom_vec.h"
#include <string.h>
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace Eigen;
using namespace std;
#define DELTA 16384
#define EPSILON 1.0e-6
enum {
LAYOUT_UNIFORM, LAYOUT_NONUNIFORM, LAYOUT_TILED
};
// several files
/* ---------------------------------------------------------------------- */
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 != 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];
sprintf(str, "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 str[128];
char line[256];
char *retpointer;
char **values;
int nwords;
// read STL solid name
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading number of triangle pairs");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords < 1) {
sprintf(str, "first line of file is incorrect");
error->one(FLERR, str);
}
// 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) {
sprintf(str, "error reading outer loop");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 2) {
sprintf(str, "error reading outer loop");
error->one(FLERR, str);
}
// read vertex lines
for (int k = 0; k < 3; k++) {
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading vertex line");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 4) {
sprintf(str, "error reading vertex line");
error->one(FLERR, str);
}
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) {
sprintf(str, "error reading endloop");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 1) {
sprintf(str, "error reading endloop");
error->one(FLERR, str);
}
// read end facet line
retpointer = fgets(line, sizeof(line), fp);
if (retpointer == NULL) {
sprintf(str, "error reading endfacet");
error->one(FLERR, str);
}
nwords = count_words(line);
if (nwords != 1) {
sprintf(str, "error reading endfacet");
error->one(FLERR, str);
}
// 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);
}
Event Timeline
Log In to Comment