Diffusion lammps (lammm-master-deprecated) Browse

fix_smd_move_triangulated_surface.cpp
354e20f43157lammm-master-deprecated

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++];
			}
			}

lammps/src/USER-SMD/fix_smd_move_triangulated_surface.cpp354e20f43157lammm-master-deprecated

fix_smd_move_triangulated_surface.cppView Options

lammps/src/USER-SMD/fix_smd_move_triangulated_surface.cpp
354e20f43157lammm-master-deprecated

fix_smd_move_triangulated_surface.cpp
View Options