Commit 892fb116 authored by sjplimp's avatar sjplimp
Browse files 

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@1608 f3b2605a-c512-4ea7-a41b-209d697bcdaa
parent 62590294

src/GRANULAR/atom_vec_granular.cpp

+1 −28
Original line number Diff line number Diff line
@@ -41,7 +41,7 @@ AtomVecGranular::AtomVecGranular(LAMMPS *lmp, int narg, char **arg) : 
  xcol_data = 5;



  atom->radius_flag = atom->density_flag = atom->rmass_flag = 1;

  atom->xorient_flag = atom->omega_flag = atom->torque_flag = 1;

  atom->omega_flag = atom->torque_flag = 1;



  PI = 4.0*atan(1.0);

}
@@ -77,8 +77,6 @@ void AtomVecGranular::grow(int n) 
  rmass = atom->rmass = (double *)

    memory->srealloc(atom->rmass,nmax*sizeof(double),"atom:rmass");



  xorient = atom->xorient = 

    memory->grow_2d_double_array(atom->xorient,nmax,3,"atom:xorient");

  omega = atom->omega = 

    memory->grow_2d_double_array(atom->omega,nmax,3,"atom:omega");

  torque = atom->torque =
@@ -107,9 +105,6 @@ void AtomVecGranular::copy(int i, int j) 
  radius[j] = radius[i];

  density[j] = density[i];

  rmass[j] = rmass[i];

  xorient[j][0] = xorient[i][0];

  xorient[j][1] = xorient[i][1];

  xorient[j][2] = xorient[i][2];

  omega[j][0] = omega[i][0];

  omega[j][1] = omega[i][1];

  omega[j][2] = omega[i][2];
@@ -408,9 +403,6 @@ int AtomVecGranular::pack_exchange(int i, double *buf) 
  buf[m++] = radius[i];

  buf[m++] = density[i];

  buf[m++] = rmass[i];

  buf[m++] = xorient[i][0];

  buf[m++] = xorient[i][1];

  buf[m++] = xorient[i][2];

  buf[m++] = omega[i][0];

  buf[m++] = omega[i][1];

  buf[m++] = omega[i][2];
@@ -445,9 +437,6 @@ int AtomVecGranular::unpack_exchange(double *buf) 
  radius[nlocal] = buf[m++];

  density[nlocal] = buf[m++];

  rmass[nlocal] = buf[m++];

  xorient[nlocal][0] = buf[m++];

  xorient[nlocal][1] = buf[m++];

  xorient[nlocal][2] = buf[m++];

  omega[nlocal][0] = buf[m++];

  omega[nlocal][1] = buf[m++];

  omega[nlocal][2] = buf[m++];
@@ -503,9 +492,6 @@ int AtomVecGranular::pack_restart(int i, double *buf) 


  buf[m++] = radius[i];

  buf[m++] = density[i];

  buf[m++] = xorient[i][0];

  buf[m++] = xorient[i][1];

  buf[m++] = xorient[i][2];

  buf[m++] = omega[i][0];

  buf[m++] = omega[i][1];

  buf[m++] = omega[i][2];
@@ -553,9 +539,6 @@ int AtomVecGranular::unpack_restart(double *buf) 
  else

    rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal];



  xorient[nlocal][0] = buf[m++];

  xorient[nlocal][1] = buf[m++];

  xorient[nlocal][2] = buf[m++];

  omega[nlocal][0] = buf[m++];

  omega[nlocal][1] = buf[m++];

  omega[nlocal][2] = buf[m++];
@@ -598,9 +581,6 @@ void AtomVecGranular::create_atom(int itype, double *coord) 
      radius[nlocal]*radius[nlocal]*radius[nlocal] * density[nlocal];

  else

    rmass[nlocal] = PI * radius[nlocal]*radius[nlocal] * density[nlocal];

  xorient[nlocal][0] = 0.0;

  xorient[nlocal][1] = 0.0;

  xorient[nlocal][2] = 0.0;

  omega[nlocal][0] = 0.0;

  omega[nlocal][1] = 0.0;

  omega[nlocal][2] = 0.0;
@@ -644,9 +624,6 @@ void AtomVecGranular::data_atom(double *coord, int imagetmp, char **values) 
  v[nlocal][0] = 0.0;

  v[nlocal][1] = 0.0;

  v[nlocal][2] = 0.0;

  xorient[nlocal][0] = 1.0;

  xorient[nlocal][1] = 0.0;

  xorient[nlocal][2] = 0.0;

  omega[nlocal][0] = 0.0;

  omega[nlocal][1] = 0.0;

  omega[nlocal][2] = 0.0;
@@ -672,9 +649,6 @@ int AtomVecGranular::data_atom_hybrid(int nlocal, char **values) 
  v[nlocal][0] = 0.0;

  v[nlocal][1] = 0.0;

  v[nlocal][2] = 0.0;

  xorient[nlocal][0] = 1.0;

  xorient[nlocal][1] = 0.0;

  xorient[nlocal][2] = 0.0;

  omega[nlocal][0] = 0.0;

  omega[nlocal][1] = 0.0;

  omega[nlocal][2] = 0.0;
@@ -727,7 +701,6 @@ double AtomVecGranular::memory_usage() 
  if (atom->memcheck("radius")) bytes += nmax * sizeof(double);

  if (atom->memcheck("density")) bytes += nmax * sizeof(double);

  if (atom->memcheck("rmass")) bytes += nmax * sizeof(double);

  if (atom->memcheck("xorient")) bytes += nmax*3 * sizeof(double);

  if (atom->memcheck("omega")) bytes += nmax*3 * sizeof(double);

  if (atom->memcheck("torque")) bytes += nmax*3 * sizeof(double);

src/GRANULAR/atom_vec_granular.h

+1 −1
Original line number Diff line number Diff line
@@ -53,7 +53,7 @@ class AtomVecGranular : public AtomVec { 
  int *tag,*type,*mask,*image;

  double **x,**v,**f;

  double *radius,*density,*rmass;

  double **xorient,**omega,**torque;

  double **omega,**torque;

};



}

src/atom.cpp

+2 −3
Original line number Diff line number Diff line
@@ -63,7 +63,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) 
  molecule = NULL;

  q = NULL;

  mu = NULL;

  xorient = quat = omega = angmom = torque = NULL;

  quat = omega = angmom = torque = NULL;

  radius = density = rmass = vfrac = NULL;



  maxspecial = 1;
@@ -89,7 +89,7 @@ Atom::Atom(LAMMPS *lmp) : Pointers(lmp) 


  molecule_flag = 0;

  q_flag = mu_flag = 0;

  xorient_flag = quat_flag = omega_flag = angmom_flag = torque_flag = 0;

  quat_flag = omega_flag = angmom_flag = torque_flag = 0;

  radius_flag = density_flag = rmass_flag = vfrac_flag = 0;



  // ntype-length arrays
@@ -153,7 +153,6 @@ Atom::~Atom() 


  memory->sfree(q);

  memory->destroy_2d_double_array(mu);

  memory->destroy_2d_double_array(xorient);

  memory->destroy_2d_double_array(quat);

  memory->destroy_2d_double_array(omega);

  memory->destroy_2d_double_array(angmom);

src/atom.h

+2 −2
Original line number Diff line number Diff line
@@ -47,7 +47,7 @@ class Atom : protected Pointers { 


  int *molecule;

  double *q,**mu;

  double **xorient,**quat,**omega,**angmom,**torque;

  double **quat,**omega,**angmom,**torque;

  double *radius,*density,*rmass,*vfrac;



  int maxspecial;
@@ -74,7 +74,7 @@ class Atom : protected Pointers { 


  int molecule_flag;

  int q_flag,mu_flag;

  int xorient_flag,quat_flag,omega_flag,angmom_flag,torque_flag;

  int quat_flag,omega_flag,angmom_flag,torque_flag;

  int radius_flag,density_flag,rmass_flag,vfrac_flag;



  // extra peratom info in restart file destined for fix & diag

src/compute_rotate_dipole.cpp→src/compute_erotate_sphere.cpp

+105 −0
Original line number Diff line number Diff line
@@ -12,8 +12,9 @@ 
------------------------------------------------------------------------- */



#include "mpi.h"

#include "compute_rotate_dipole.h"

#include "compute_erotate_sphere.h"

#include "atom.h"

#include "force.h"

#include "domain.h"

#include "group.h"

#include "error.h"
@@ -21,70 +22,84 @@ 
using namespace LAMMPS_NS;



#define INVOKED_SCALAR 1

#define INERTIA3D 0.4          // moments of inertia for sphere and disk

#define INERTIA2D 0.5



#define INERTIA 0.4          // moment of inertia for sphere



/* ---------------------------------------------------------------------- */



ComputeRotateDipole::ComputeRotateDipole(LAMMPS *lmp, int narg, char **arg) : 

ComputeERotateSphere::ComputeERotateSphere(LAMMPS *lmp, int narg, char **arg) :

  Compute(lmp, narg, arg)

{

  if (narg != 3) error->all("Illegal compute rotate/dipole command");

  if (narg != 3) error->all("Illegal compute erotate/sphere command");



  if (atom->dipole == NULL || !atom->omega_flag)

    error->all("Compute rotate/dipole requires atom attributes dipole, omega");

  if (!atom->omega_flag) 

    error->all("Compute erotate/sphere requires atom attribute omega");



  scalar_flag = 1;

  extscalar = 1;



  inertia = NULL;

  inertia = new double[atom->ntypes+1];

}



/* ---------------------------------------------------------------------- */



ComputeRotateDipole::~ComputeRotateDipole()

ComputeERotateSphere::~ComputeERotateSphere()

{

  delete [] inertia;

}



/* ---------------------------------------------------------------------- */



void ComputeRotateDipole::init()

void ComputeERotateSphere::init()

{

  delete [] inertia;

  inertia = new double[atom->ntypes+1];

  pfactor = 0.5 * force->mvv2e * INERTIA;



  if (atom->mass && !atom->shape)

    error->all("Compute erotate/sphere requires atom attribute shape");

  if (!atom->mass && (!atom->radius_flag || !atom->rmass_flag))

    error->all("Compute erotate/sphere requires atom attributes radius, rmass");



  if (atom->mass) {

    double *mass = atom->mass;

    double **shape = atom->shape;

    

  if (domain->dimension == 3)

    for (int i = 1; i <= atom->ntypes; i++)

      inertia[i] = INERTIA3D * mass[i] * 0.25*shape[i][0]*shape[i][0];

  else

    for (int i = 1; i <= atom->ntypes; i++)

      inertia[i] = INERTIA2D * mass[i] * 0.25*shape[i][0]*shape[i][0];

    for (int i = 1; i <= atom->ntypes; i++) {

      if (shape[i][0] != shape[i][1] || shape[i][0] != shape[i][2])

	error->all("Compute rotate requires spherical particle shapes");

      inertia[i] = 0.25*shape[i][0]*shape[i][0] * mass[i];

    }

  }

}



/* ---------------------------------------------------------------------- */



double ComputeRotateDipole::compute_scalar()

double ComputeERotateSphere::compute_scalar()

{

  invoked |= INVOKED_SCALAR;



  double *dipole = atom->dipole;

  double **omega = atom->omega;

  int *type = atom->type;

  double *radius = atom->radius;

  double *rmass = atom->rmass;

  double *mass = atom->mass;

  int *mask = atom->mask;

  int *type = atom->type;

  int nlocal = atom->nlocal;



  double erot = 0.0;

  double erotate = 0.0;



  if (mass) {

    for (int i = 0; i < nlocal; i++)

      if (mask[i] & groupbit)

	erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 

		    omega[i][2]*omega[i][2]) * inertia[type[i]];

  } else {

    for (int i = 0; i < nlocal; i++) 

    if (mask[i] & groupbit && dipole[type[i]] > 0.0)

      erot += inertia[type[i]] * 

	(omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 

	 omega[i][2]*omega[i][2]);

      if (mask[i] & groupbit)

	erotate += (omega[i][0]*omega[i][0] + omega[i][1]*omega[i][1] + 

		    omega[i][2]*omega[i][2]) * radius[i]*radius[i]*rmass[i];

  }



  MPI_Allreduce(&erot,&scalar,1,MPI_DOUBLE,MPI_SUM,world);

  scalar *= 0.5;

  MPI_Allreduce(&erotate,&scalar,1,MPI_DOUBLE,MPI_SUM,world);

  scalar *= pfactor;

  return scalar;

}

CRA Git | Maintained and supported by SUSTech CRA and CCSE