Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F64630572
fix_smd_move_triangulated_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
Tue, May 28, 07:17
Size
13 KB
Mime Type
text/x-c
Expires
Thu, May 30, 07:17 (2 d)
Engine
blob
Format
Raw Data
Handle
17790605
Attached To
rLAMMPS lammps
fix_smd_move_triangulated_surface.cpp
View Options
/* ----------------------------------------------------------------------
*
* *** Smooth Mach Dynamics ***
*
* This file is part of the USER-SMD package for LAMMPS.
* Copyright (2014) Georg C. Ganzenmueller, georg.ganzenmueller@emi.fhg.de
* Fraunhofer Ernst-Mach Institute for High-Speed Dynamics, EMI,
* Eckerstrasse 4, D-79104 Freiburg i.Br, Germany.
*
* ----------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
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.
------------------------------------------------------------------------- */
#include <stdio.h>
#include <string.h>
#include "fix_smd_move_triangulated_surface.h"
#include <math.h>
#include <stdlib.h>
#include <string.h>
#include "atom.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "update.h"
#include "integrate.h"
#include "respa.h"
#include "memory.h"
#include "error.h"
#include "pair.h"
#include "domain.h"
#include <Eigen/Eigen>
#include "math_const.h"
using namespace Eigen;
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
using namespace std;
/* ---------------------------------------------------------------------- */
FixSMDMoveTriSurf::FixSMDMoveTriSurf(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg) {
if (atom->smd_flag != 1) {
error->all(FLERR, "fix fix smd/move_tri_surf command requires atom_style smd");
}
if (narg < 3)
error->all(FLERR, "Illegal number of arguments for fix fix smd/move_tri_surf command");
rotateFlag = linearFlag = wiggleFlag = false;
wiggle_direction = 1.0;
wiggle_max_travel = 0.0;
int iarg = 3;
if (comm->me == 0) {
printf("\n>>========>>========>>========>>========>>========>>========>>========>>========\n");
printf("fix fix smd/move_tri_surf is active for group: %s \n", arg[1]);
}
while (true) {
if (iarg >= narg) {
break;
}
if (strcmp(arg[iarg], "*LINEAR") == 0) {
linearFlag = true;
if (comm->me == 0) {
printf("... will move surface in a linear fashion\n");
}
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three floats for velocity following *LINEAR");
}
vx = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three floats for velocity following *LINEAR");
}
vy = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected three floats for velocity following *LINEAR");
}
vz = force->numeric(FLERR, arg[iarg]);
} else if (strcmp(arg[iarg], "*WIGGLE") == 0) {
wiggleFlag = true;
if (comm->me == 0) {
printf("... will move surface in wiggle fashion\n");
}
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
vx = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
vy = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
vz = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 4 floats following *WIGGLE : vx vy vz max_travel");
}
wiggle_max_travel = force->numeric(FLERR, arg[iarg]);
} else if (strcmp(arg[iarg], "*ROTATE") == 0) {
rotateFlag = true;
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
origin(0) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
origin(1) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
origin(2) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_axis(0) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_axis(1) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_axis(2) = force->numeric(FLERR, arg[iarg]);
iarg++;
if (iarg == narg) {
error->all(FLERR, "expected 7 floats following *ROTATE: origin, rotation axis, and rotation period");
}
rotation_period = force->numeric(FLERR, arg[iarg]);
/*
* construct rotation matrix
*/
u_cross(0, 0) = 0.0;
u_cross(0, 1) = -rotation_axis(2);
u_cross(0, 2) = rotation_axis(1);
u_cross(1, 0) = rotation_axis(2);
u_cross(1, 1) = 0.0;
u_cross(1, 2) = -rotation_axis(0);
u_cross(2, 0) = -rotation_axis(1);
u_cross(2, 1) = rotation_axis(0);
u_cross(2, 2) = 0.0;
uxu = rotation_axis * rotation_axis.transpose();
if (comm->me == 0) {
printf("will rotate with period %f\n", rotation_period);
}
} else {
char msg[128];
sprintf(msg, "Illegal keyword for fix smd/move_tri_surf: %s\n", arg[iarg]);
error->all(FLERR, msg);
}
iarg++;
}
if (comm->me == 0) {
printf(">>========>>========>>========>>========>>========>>========>>========>>========\n\n");
}
// set comm sizes needed by this fix
comm_forward = 12;
//atom->add_callback(0);
//atom->add_callback(1);
time_integrate = 1;
}
/* ---------------------------------------------------------------------- */
FixSMDMoveTriSurf::~FixSMDMoveTriSurf()
{
// unregister callbacks to this fix from Atom class
//atom->delete_callback(id,0);
//atom->delete_callback(id,1);
}
/* ---------------------------------------------------------------------- */
int FixSMDMoveTriSurf::setmask() {
int mask = 0;
mask |= INITIAL_INTEGRATE;
//mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixSMDMoveTriSurf::init() {
dtv = update->dt;
}
/* ----------------------------------------------------------------------
------------------------------------------------------------------------- */
void FixSMDMoveTriSurf::initial_integrate(int vflag) {
double **x = atom->x;
double **x0 = atom->x0;
double **v = atom->v;
double **vest = atom->vest;
double **smd_data_9 = atom->smd_data_9;
tagint *mol = atom->molecule;
int *mask = atom->mask;
int nlocal = atom->nlocal;
double phi;
int i;
Matrix3d eye, Rot;
eye.setIdentity();
Vector3d v1, v2, v3, n, point, rotated_point, R, vel;
if (igroup == atom->firstgroup)
nlocal = atom->nfirst;
if (linearFlag) { // translate particles
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] = vx;
v[i][1] = vy;
v[i][2] = vz;
vest[i][0] = vx;
vest[i][1] = vy;
vest[i][2] = vz;
x[i][0] += dtv * vx;
x[i][1] += dtv * vy;
x[i][2] += dtv * vz;
/*
* if this is a triangle, move the vertices as well
*/
if (mol[i] >= 65535) {
smd_data_9[i][0] += dtv * vx;
smd_data_9[i][1] += dtv * vy;
smd_data_9[i][2] += dtv * vz;
smd_data_9[i][3] += dtv * vx;
smd_data_9[i][4] += dtv * vy;
smd_data_9[i][5] += dtv * vz;
smd_data_9[i][6] += dtv * vx;
smd_data_9[i][7] += dtv * vy;
smd_data_9[i][8] += dtv * vz;
}
}
}
}
if (wiggleFlag) { // wiggle particles forward and backward
wiggle_travel += sqrt(vx * vx + vy * vy + vz * vz ) * dtv;
double wiggle_vx = wiggle_direction * vx;
double wiggle_vy = wiggle_direction * vy;
double wiggle_vz = wiggle_direction * vz;
//printf("wiggle vz is %f, wiggle_max_travel is %f, dir=%f\n", wiggle_vz, wiggle_max_travel, wiggle_direction);
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
v[i][0] = wiggle_vx;
v[i][1] = wiggle_vy;
v[i][2] = wiggle_vz;
vest[i][0] = wiggle_vx;
vest[i][1] = wiggle_vy;
vest[i][2] = wiggle_vz;
x[i][0] += dtv * wiggle_vx;
x[i][1] += dtv * wiggle_vy;
x[i][2] += dtv * wiggle_vz;
/*
* if this is a triangle, move the vertices as well
*/
if (mol[i] >= 65535) {
smd_data_9[i][0] += dtv * wiggle_vx;
smd_data_9[i][1] += dtv * wiggle_vy;
smd_data_9[i][2] += dtv * wiggle_vz;
smd_data_9[i][3] += dtv * wiggle_vx;
smd_data_9[i][4] += dtv * wiggle_vy;
smd_data_9[i][5] += dtv * wiggle_vz;
smd_data_9[i][6] += dtv * wiggle_vx;
smd_data_9[i][7] += dtv * wiggle_vy;
smd_data_9[i][8] += dtv * wiggle_vz;
}
}
}
if (wiggle_travel >= wiggle_max_travel) {
wiggle_direction *= -1.0;
wiggle_travel = 0.0;
}
}
if (rotateFlag) { // rotate particles
Vector3d xnew, R_new, x_correct;
/*
* rotation angle and matrix form of Rodrigues' rotation formula
*/
phi = MY_2PI * dtv / rotation_period;
//printf("dt=%f, phi =%f, T=%f\n", dtv, phi, rotation_period);
Rot = cos(phi) * eye + sin(phi) * u_cross + (1.0 - cos(phi)) * uxu;
for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
/*
* generate vector R from origin to point which is to be rotated
*/
point << x[i][0], x[i][1], x[i][2];
R = point - origin;
/*
* rotate vector and shift away from origin
*/
rotated_point = Rot * R + origin;
/*
* determine velocity
*/
vel = (rotated_point - point) / update->dt;
/*
* assign new velocities and coordinates
*/
v[i][0] = vel(0);
v[i][1] = vel(1);
v[i][2] = vel(2);
vest[i][0] = vel(0);
vest[i][1] = vel(1);
vest[i][2] = vel(2);
x[i][0] = rotated_point(0);
x[i][1] = rotated_point(1);
x[i][2] = rotated_point(2);
/*
* if this is a triangle, rotate the vertices as well
*/
if (mol[i] >= 65535) {
v1 << smd_data_9[i][0], smd_data_9[i][1], smd_data_9[i][2];
R = v1 - origin;
rotated_point = Rot * R + origin;
vel = (rotated_point - v1) / update->dt;
smd_data_9[i][0] = rotated_point(0);
smd_data_9[i][1] = rotated_point(1);
smd_data_9[i][2] = rotated_point(2);
v1 = rotated_point;
v2 << smd_data_9[i][3], smd_data_9[i][4], smd_data_9[i][5];
R = v2 - origin;
rotated_point = Rot * R + origin;
vel = (rotated_point - v2) / update->dt;
smd_data_9[i][3] = rotated_point(0);
smd_data_9[i][4] = rotated_point(1);
smd_data_9[i][5] = rotated_point(2);
v2 = rotated_point;
v3 << smd_data_9[i][6], smd_data_9[i][7], smd_data_9[i][8];
R = v3 - origin;
rotated_point = Rot * R + origin;
vel = (rotated_point - v3) / update->dt;
smd_data_9[i][6] = rotated_point(0);
smd_data_9[i][7] = rotated_point(1);
smd_data_9[i][8] = rotated_point(2);
v3 = rotated_point;
// recalculate triangle normal
n = (v2 - v1).cross(v2 - v3);
n /= n.norm();
x0[i][0] = n(0);
x0[i][1] = n(1);
x0[i][2] = n(2);
}
}
}
}
// we changed smd_data_9, x0. perform communication to ghosts
comm->forward_comm_fix(this);
}
/* ---------------------------------------------------------------------- */
void FixSMDMoveTriSurf::reset_dt() {
dtv = update->dt;
}
/* ---------------------------------------------------------------------- */
int FixSMDMoveTriSurf::pack_forward_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) {
int i, j, m;
double **x0 = atom->x0;
double **smd_data_9 = atom->smd_data_9;
//printf("in FixSMDIntegrateTlsph::pack_forward_comm\n");
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
buf[m++] = x0[j][0];
buf[m++] = x0[j][1];
buf[m++] = x0[j][2];
buf[m++] = smd_data_9[j][0];
buf[m++] = smd_data_9[j][1];
buf[m++] = smd_data_9[j][2];
buf[m++] = smd_data_9[j][3];
buf[m++] = smd_data_9[j][4];
buf[m++] = smd_data_9[j][5];
buf[m++] = smd_data_9[j][6];
buf[m++] = smd_data_9[j][7];
buf[m++] = smd_data_9[j][8];
}
return m;
}
/* ---------------------------------------------------------------------- */
void FixSMDMoveTriSurf::unpack_forward_comm(int n, int first, double *buf) {
int i, m, last;
double **x0 = atom->x0;
double **smd_data_9 = atom->smd_data_9;
//printf("in FixSMDMoveTriSurf::unpack_forward_comm\n");
m = 0;
last = first + n;
for (i = first; i < last; i++) {
x0[i][0] = buf[m++];
x0[i][1] = buf[m++];
x0[i][2] = buf[m++];
smd_data_9[i][0] = buf[m++];
smd_data_9[i][1] = buf[m++];
smd_data_9[i][2] = buf[m++];
smd_data_9[i][3] = buf[m++];
smd_data_9[i][4] = buf[m++];
smd_data_9[i][5] = buf[m++];
smd_data_9[i][6] = buf[m++];
smd_data_9[i][7] = buf[m++];
smd_data_9[i][8] = buf[m++];
}
}
Event Timeline
Log In to Comment