diff --git a/src/PERI/compute_damage_atom.cpp b/src/PERI/compute_damage_atom.cpp index 9b51cd34b3..7a8cbfd3dd 100644 --- a/src/PERI/compute_damage_atom.cpp +++ b/src/PERI/compute_damage_atom.cpp @@ -100,23 +100,22 @@ void ComputeDamageAtom::compute_peratom() if (mask[i] & groupbit) { jnum = npartner[i]; damage_temp = 0.0; - for (jj = 0; jj < jnum; jj++) { if (partner[i][jj] == 0) continue; - + // look up local index of this partner particle // skip if particle is "lost" - + j = atom->map(partner[i][jj]); if (j < 0) continue; - + damage_temp += vfrac[j]; } - } - else damage_temp = vinter[i]; - if (vinter[i] != 0.0) damage[i] = 1.0 - damage_temp/vinter[i]; - else damage[i] = 0.0; + if (vinter[i] != 0.0) damage[i] = 1.0 - damage_temp/vinter[i]; + else damage[i] = 0.0; + + } else damage[i] = 0.0; } } diff --git a/src/PERI/compute_dilatation_atom.cpp b/src/PERI/compute_dilatation_atom.cpp new file mode 100644 index 0000000000..366c5bcda9 --- /dev/null +++ b/src/PERI/compute_dilatation_atom.cpp @@ -0,0 +1,127 @@ +/* ---------------------------------------------------------------------- + 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 author: Rezwanur Rahman, John Foster (UTSA) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_dilatation_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "comm.h" +#include "force.h" +#include "pair.h" +#include "pair_peri_lps.h" +#include "pair_peri_pmb.h" +#include "pair_peri_ves.h" +#include "pair_peri_eps.h" +#include "fix_peri_neigh.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputeDilatationAtom:: +ComputeDilatationAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all(FLERR,"Illegal compute Dilatation/atom command"); + + peratom_flag = 1; + size_peratom_cols = 0; + + nmax = 0; + dilatation = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputeDilatationAtom::~ComputeDilatationAtom() +{ + memory->destroy(dilatation); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDilatationAtom::init() +{ + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"dilatation/peri") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one compute dilatation/atom"); + + // check PD pair style + + isPMB = isLPS = isVES = isEPS = 0; + if (force->pair_match("peri/pmb",1)) isPMB = 1; + if (force->pair_match("peri/lps",1)) isLPS = 1; + if (force->pair_match("peri/ves",1)) isVES = 1; + if (force->pair_match("peri/eps",1)) isEPS = 1; + + if (isPMB) + error->all(FLERR,"Compute dilatation/atom cannot be used " + "with this pair style"); + + // find associated PERI_NEIGH fix that must exist + + int ifix_peri = -1; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + if (ifix_peri == -1) + error->all(FLERR,"Compute dilatation/atom requires Peridynamic pair style"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeDilatationAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow dilatation array if necessary + + if (atom->nlocal > nmax) { + memory->destroy(dilatation); + nmax = atom->nmax; + memory->create(dilatation,nmax,"dilatation/atom:dilatation"); + vector_atom = dilatation; + } + + // extract dilatation for each atom in group + + double *theta; + Pair *anypair = force->pair_match("peri",0); + if (isLPS) theta = ((PairPeriLPS *) anypair)->theta; + if (isVES) theta = ((PairPeriVES *) anypair)->theta; + if (isEPS) theta = ((PairPeriEPS *) anypair)->theta; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) dilatation[i] = theta[i]; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputeDilatationAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/PERI/compute_dilatation_atom.h b/src/PERI/compute_dilatation_atom.h new file mode 100644 index 0000000000..5e839cc3d0 --- /dev/null +++ b/src/PERI/compute_dilatation_atom.h @@ -0,0 +1,48 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(dilatation/atom,ComputeDilatationAtom) + +#else + +#ifndef LMP_COMPUTE_DILATATION_ATOM_H +#define LMP_COMPUTE_DILATATION_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeDilatationAtom : public Compute { + friend class PairPeriPMB; + friend class PairPeriLPS; + friend class PairPeriVES; + friend class PairPeriEPS; + public: + ComputeDilatationAtom(class LAMMPS *, int, char **); + ~ComputeDilatationAtom(); + void init(); + void compute_peratom(); + double memory_usage(); + + private: + int nmax; + double *dilatation; + int isPMB,isLPS,isVES,isEPS; +}; + +} + +#endif +#endif diff --git a/src/PERI/compute_plasticity_atom.cpp b/src/PERI/compute_plasticity_atom.cpp new file mode 100644 index 0000000000..ac089fe7a6 --- /dev/null +++ b/src/PERI/compute_plasticity_atom.cpp @@ -0,0 +1,111 @@ +/* ---------------------------------------------------------------------- + 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 author: Rezwanur Rahman, John Foster (UTSA) +------------------------------------------------------------------------- */ + +#include "string.h" +#include "compute_plasticity_atom.h" +#include "atom.h" +#include "update.h" +#include "modify.h" +#include "comm.h" +#include "force.h" +#include "pair_peri_pmb.h" +#include "fix_peri_neigh.h" +#include "memory.h" +#include "error.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +ComputePlasticityAtom:: +ComputePlasticityAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg) +{ + if (narg != 3) error->all(FLERR,"Illegal compute plasticity/atom command"); + + if (!force->pair_match("peri/eps",1)) + error->all(FLERR,"Compute plasticity/atom cannot be used " + "with this pair style"); + + peratom_flag = 1; + size_peratom_cols = 0; + + nmax = 0; + plasticity = NULL; +} + +/* ---------------------------------------------------------------------- */ + +ComputePlasticityAtom::~ComputePlasticityAtom() +{ + memory->destroy(plasticity); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePlasticityAtom::init() +{ + int count = 0; + for (int i = 0; i < modify->ncompute; i++) + if (strcmp(modify->compute[i]->style,"plasticity/peri") == 0) count++; + if (count > 1 && comm->me == 0) + error->warning(FLERR,"More than one compute plasticity/atom"); + + // find associated PERI_NEIGH fix that must exist + + ifix_peri = -1; + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + if (ifix_peri == -1) + error->all(FLERR,"Compute plasticity/atom requires Peridynamic pair style"); +} + +/* ---------------------------------------------------------------------- */ + +void ComputePlasticityAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + // grow damage array if necessary + + if (atom->nlocal > nmax) { + memory->destroy(plasticity); + nmax = atom->nmax; + memory->create(plasticity,nmax,"plasticity/atom:plasticity"); + vector_atom = plasticity; + } + + // extract plasticity for each atom in group + + double *lambdaValue = ((FixPeriNeigh *) modify->fix[ifix_peri])->lambdaValue; + + int *mask = atom->mask; + int nlocal = atom->nlocal; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) plasticity[i] = lambdaValue[i]; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based array +------------------------------------------------------------------------- */ + +double ComputePlasticityAtom::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} diff --git a/src/PERI/compute_plasticity_atom.h b/src/PERI/compute_plasticity_atom.h new file mode 100644 index 0000000000..9800b1253b --- /dev/null +++ b/src/PERI/compute_plasticity_atom.h @@ -0,0 +1,44 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS + +ComputeStyle(plasticity/atom,ComputePlasticityAtom) + +#else + +#ifndef LMP_COMPUTE_PLASTICITY_ATOM_H +#define LMP_COMPUTE_PLASTICITY_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputePlasticityAtom : public Compute { + public: + ComputePlasticityAtom(class LAMMPS *, int, char **); + ~ComputePlasticityAtom(); + void init(); + void compute_peratom(); + double memory_usage(); + + private: + int nmax; + double *plasticity; + int ifix_peri; +}; + +} + +#endif +#endif diff --git a/src/PERI/fix_peri_neigh.cpp b/src/PERI/fix_peri_neigh.cpp index 23e0a392af..190bc0460d 100644 --- a/src/PERI/fix_peri_neigh.cpp +++ b/src/PERI/fix_peri_neigh.cpp @@ -20,6 +20,7 @@ #include "pair_peri_pmb.h" #include "pair_peri_lps.h" #include "pair_peri_ves.h" +#include "pair_peri_eps.h" #include "atom.h" #include "domain.h" #include "force.h" @@ -41,10 +42,11 @@ using namespace FixConst; FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : Fix(lmp, narg, arg) { - isPMB = isLPS = isVES = 0; + isPMB = isLPS = isVES = isEPS = 0; if (force->pair_match("peri/pmb",1)) isPMB = 1; if (force->pair_match("peri/lps",1)) isLPS = 1; if (force->pair_match("peri/ves",1)) isVES = 1; + if (force->pair_match("peri/eps",1)) isEPS = 1; restart_global = 1; restart_peratom = 1; @@ -59,6 +61,8 @@ FixPeriNeigh::FixPeriNeigh(LAMMPS *lmp,int narg, char **arg) : partner = NULL; deviatorextention = NULL; deviatorBackextention = NULL; + deviatorPlasticextension = NULL; + lambdaValue = NULL; r0 = NULL; vinter = NULL; wvolume = NULL; @@ -92,6 +96,8 @@ FixPeriNeigh::~FixPeriNeigh() memory->destroy(partner); memory->destroy(deviatorextention); memory->destroy(deviatorBackextention); + memory->destroy(deviatorPlasticextension); + memory->destroy(lambdaValue); memory->destroy(r0); memory->destroy(vinter); memory->destroy(wvolume); @@ -218,6 +224,8 @@ void FixPeriNeigh::setup(int vflag) memory->destroy(partner); memory->destroy(deviatorextention); memory->destroy(deviatorBackextention); + memory->destroy(deviatorPlasticextension); + memory->destroy(lambdaValue); memory->destroy(r0); memory->destroy(npartner); @@ -225,6 +233,8 @@ void FixPeriNeigh::setup(int vflag) partner = NULL; deviatorextention = NULL; deviatorBackextention = NULL; + deviatorPlasticextension = NULL; + lambdaValue = NULL; r0 = NULL; grow_arrays(atom->nmax); @@ -235,6 +245,7 @@ void FixPeriNeigh::setup(int vflag) npartner[i] = 0; vinter[i] = 0.0; wvolume[i] = 0.0; + if (isEPS) lambdaValue[i] = 0.0; } for (ii = 0; ii < inum; ii++) { @@ -261,6 +272,8 @@ void FixPeriNeigh::setup(int vflag) if (isVES) deviatorextention[i][npartner[i]] = deviatorBackextention[i][npartner[i]] = 0.0; + if (isEPS) + deviatorPlasticextension[i][npartner[i]] = 0.0; r0[i][npartner[i]] = sqrt(rsq); npartner[i]++; vinter[i] += vfrac[j]; @@ -292,6 +305,7 @@ void FixPeriNeigh::setup(int vflag) PairPeriLPS *pairlps = static_cast(anypair); PairPeriPMB *pairpmb = static_cast(anypair); PairPeriVES *pairves = static_cast(anypair); + PairPeriEPS *paireps = static_cast(anypair); for (i = 0; i < nlocal; i++) { double xtmp0 = x0[i][0]; @@ -319,6 +333,7 @@ void FixPeriNeigh::setup(int vflag) double delx0 = xtmp0 - x0[j][0]; double dely0 = ytmp0 - x0[j][1]; double delz0 = ztmp0 - x0[j][2]; + double rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; jtype = type[j]; @@ -332,7 +347,6 @@ void FixPeriNeigh::setup(int vflag) else vfrac_scale = 1.0; // for PMB, influence = 1.0, otherwise invoke influence function - if (isPMB) wvolume[i] += 1.0 * rsq0 * vfrac[j] * vfrac_scale; else if (isLPS) @@ -341,6 +355,9 @@ void FixPeriNeigh::setup(int vflag) else if (isVES) wvolume[i] += pairves->influence_function(delx0,dely0,delz0) * rsq0 * vfrac[j] * vfrac_scale; + else if (isEPS) + wvolume[i] += paireps->influence_function(delx0,dely0,delz0) * + rsq0 * vfrac[j] * vfrac_scale; } } @@ -377,12 +394,16 @@ double FixPeriNeigh::memory_usage() { int nmax = atom->nmax; int bytes = nmax * sizeof(int); - bytes += nmax*maxpartner * sizeof(int); + bytes += nmax*maxpartner * sizeof(tagint); bytes += nmax*maxpartner * sizeof(double); if (isVES) { bytes += nmax*maxpartner * sizeof(double); bytes += nmax*maxpartner * sizeof(double); } + if (isEPS) { + bytes += nmax*maxpartner * sizeof(double); + bytes += nmax * sizeof(double); + } bytes += nmax * sizeof(double); bytes += nmax * sizeof(double); return bytes; @@ -402,7 +423,10 @@ void FixPeriNeigh::grow_arrays(int nmax) memory->grow(deviatorBackextention,nmax,maxpartner, "peri_neigh:deviatorBackextention"); } + if (isEPS) memory->grow(deviatorPlasticextension,nmax,maxpartner, + "peri_neigh:deviatorPlasticextension"); memory->grow(r0,nmax,maxpartner,"peri_neigh:r0"); + if (isEPS) memory->grow(lambdaValue,nmax,"peri_neigh:lambdaValue"); memory->grow(vinter,nmax,"peri_neigh:vinter"); memory->grow(wvolume,nmax,"peri_neigh:wvolume"); } @@ -420,8 +444,11 @@ void FixPeriNeigh::copy_arrays(int i, int j, int delflag) deviatorextention[j][m] = deviatorextention[i][m]; deviatorBackextention[j][m] = deviatorBackextention[i][m]; } + if (isEPS) + deviatorPlasticextension[j][m] = deviatorPlasticextension[i][m]; r0[j][m] = r0[i][m]; } + if (isEPS) lambdaValue[j] = lambdaValue[i]; vinter[j] = vinter[i]; wvolume[j] = wvolume[i]; } @@ -443,10 +470,13 @@ int FixPeriNeigh::pack_exchange(int i, double *buf) buf[m++] = deviatorextention[i][n]; buf[m++] = deviatorBackextention[i][n]; } + if (isEPS) buf[m++] = deviatorPlasticextension[i][n]; buf[m++] = r0[i][n]; } if (isVES) buf[0] = m/4; - else buf[0] = m/2; + else if (isEPS) buf[0] = m/3; + else buf[0] = m/2; + if (isEPS) buf[m++] = lambdaValue[i]; buf[m++] = vinter[i]; buf[m++] = wvolume[i]; return m; @@ -466,8 +496,10 @@ int FixPeriNeigh::unpack_exchange(int nlocal, double *buf) deviatorextention[nlocal][n] = buf[m++]; deviatorBackextention[nlocal][n] = buf[m++]; } + if (isEPS) deviatorPlasticextension[nlocal][n] = buf[m++]; r0[nlocal][n] = buf[m++]; } + if (isEPS) lambdaValue[nlocal] = buf[m++]; vinter[nlocal] = buf[m++]; wvolume[nlocal] = buf[m++]; return m; @@ -543,8 +575,8 @@ int FixPeriNeigh::pack_restart(int i, double *buf) { int m = 0; if (isVES) buf[m++] = 4*npartner[i] + 4; + else if (isEPS) buf[m++] = 3*npartner[i] + 5; else buf[m++] = 2*npartner[i] + 4; - buf[m++] = npartner[i]; for (int n = 0; n < npartner[i]; n++) { buf[m++] = partner[i][n]; @@ -552,8 +584,10 @@ int FixPeriNeigh::pack_restart(int i, double *buf) buf[m++] = deviatorextention[i][n]; buf[m++] = deviatorBackextention[i][n]; } + if (isEPS) buf[m++] = deviatorPlasticextension[i][n]; buf[m++] = r0[i][n]; } + if (isEPS) buf[m++] = lambdaValue[i]; buf[m++] = vinter[i]; buf[m++] = wvolume[i]; return m; @@ -581,8 +615,10 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth) deviatorextention[nlocal][n] = extra[nlocal][m++]; deviatorBackextention[nlocal][n] = extra[nlocal][m++]; } + if (isEPS) deviatorPlasticextension[nlocal][n] = extra[nlocal][m++]; r0[nlocal][n] = extra[nlocal][m++]; } + if (isEPS) lambdaValue[nlocal] = extra[nlocal][m++]; vinter[nlocal] = extra[nlocal][m++]; wvolume[nlocal] = extra[nlocal][m++]; } @@ -594,6 +630,7 @@ void FixPeriNeigh::unpack_restart(int nlocal, int nth) int FixPeriNeigh::maxsize_restart() { if (isVES) return 4*maxpartner + 4; + if (isEPS) return 3*maxpartner + 5; return 2*maxpartner + 4; } @@ -604,5 +641,6 @@ int FixPeriNeigh::maxsize_restart() int FixPeriNeigh::size_restart(int nlocal) { if (isVES) return 4*npartner[nlocal] + 4; + if (isEPS) return 3*npartner[nlocal] + 5; return 2*npartner[nlocal] + 4; } diff --git a/src/PERI/fix_peri_neigh.h b/src/PERI/fix_peri_neigh.h index bfa815be6b..1c493457c6 100644 --- a/src/PERI/fix_peri_neigh.h +++ b/src/PERI/fix_peri_neigh.h @@ -28,9 +28,11 @@ class FixPeriNeigh : public Fix { friend class PairPeriPMB; friend class PairPeriPMBOMP; friend class PairPeriLPS; - friend class PairPeriVES; //NEW + friend class PairPeriVES; + friend class PairPeriEPS; friend class PairPeriLPSOMP; friend class ComputeDamageAtom; + friend class ComputePlasticityAtom; public: FixPeriNeigh(class LAMMPS *,int, char **); @@ -63,14 +65,14 @@ class FixPeriNeigh : public Fix { tagint **partner; // neighs for each atom, stored as global IDs double **deviatorextention; // Deviatoric extention double **deviatorBackextention; // Deviatoric back extention + double **deviatorPlasticextension; // Deviatoric plastic extension + double *lambdaValue; double **r0; // initial distance to partners - double **r1; // Instanteneous distance to partners *** NEW *** - double *thetaOld; // Dilatation Old one + double **r1; // instanteneous distance to partners + double *thetaValue; // dilatation double *vinter; // sum of vfrac for bonded neighbors double *wvolume; // weighted volume of particle - int isPMB; - int isLPS; - int isVES; + int isPMB,isLPS,isVES,isEPS; // which flavor of PD class NeighList *list; }; diff --git a/src/PERI/pair_peri_eps.cpp b/src/PERI/pair_peri_eps.cpp new file mode 100644 index 0000000000..fe03ac53a1 --- /dev/null +++ b/src/PERI/pair_peri_eps.cpp @@ -0,0 +1,843 @@ +/* ---------------------------------------------------------------------- + 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 author: Rezwanur Rahman, John Foster (UTSA) +------------------------------------------------------------------------- */ + +#include "math.h" +#include "stdlib.h" +#include "string.h" +#include "pair_peri_eps.h" +#include "atom.h" +#include "domain.h" +#include "lattice.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_peri_neigh.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "memory.h" +#include "error.h" +#include "update.h" + +using namespace LAMMPS_NS; + +/* ---------------------------------------------------------------------- */ + +PairPeriEPS::PairPeriEPS(LAMMPS *lmp) : Pair(lmp) +{ + for (int i = 0; i < 6; i++) virial[i] = 0.0; + no_virial_fdotr_compute = 1; + single_enable = 0; + + ifix_peri = -1; + + nmax = 0; + s0_new = NULL; + theta = NULL; + + bulkmodulus = NULL; + shearmodulus = NULL; + s00 = alpha = NULL; + cut = NULL; + m_yieldstress = NULL; + + // set comm size needed by this Pair + // comm_reverse not needed + + comm_forward = 1; +} + +/* ---------------------------------------------------------------------- */ + +PairPeriEPS::~PairPeriEPS() +{ + if (ifix_peri >= 0) modify->delete_fix("PERI_NEIGH"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(bulkmodulus); + memory->destroy(shearmodulus); + memory->destroy(s00); + memory->destroy(alpha); + memory->destroy(cut); + memory->destroy(m_yieldstress); + memory->destroy(theta); + memory->destroy(s0_new); + } +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriEPS::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0,rsq0; + double rsq,r,dr,dr1,rk,rkNew,evdwl,fpair,fbond; + double ed,fbondElastoPlastic,fbondFinal; + double deltalambda,edpNp1; + int *ilist,*jlist,*numneigh,**firstneigh; + double d_ij,delta,stretch; + + evdwl = 0.0; + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = eflag_global = eflag_atom = 0; + + double **f = atom->f; + double **x = atom->x; + int *type = atom->type; + int nlocal = atom->nlocal; + + double timestepsize = update->dt; + double *vfrac = atom->vfrac; + double *s0 = atom->s0; + double **x0 = atom->x0; + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + double **deviatorPlasticextension = + ((FixPeriNeigh *) modify->fix[ifix_peri])->deviatorPlasticextension; + tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + double *lambdaValue = ((FixPeriNeigh *) modify->fix[ifix_peri])->lambdaValue; + + // lc = lattice constant + // init_style guarantees it's the same in x, y, and z + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + double vfrac_scale = 1.0; + + // short-range forces + + int newton_pair = force->newton_pair; + int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + + // loop over neighbors of my atoms + // need minimg() for x0 difference since not ghosted + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jlist = firstneigh[i]; + jnum = numneigh[i]; + + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= NEIGHMASK; + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + rsq0 = delx0*delx0 + dely0*dely0 + delz0*delz0; + jtype = type[j]; + + r = sqrt(rsq); + + // short-range interaction distance based on initial particle position + // 0.9 and 1.35 are constants + + d_ij = MIN(0.9*sqrt(rsq0),1.35*lc); + + // short-range contact forces + // 15 is constant taken from the EMU Theory Manual + // Silling, 12 May 2005, p 18 + + if (r < d_ij) { + dr = r - d_ij; + + // kshort based upon short-range force constant + // of the bond-based theory used in PMB model + + double kshort = (15.0 * 18.0 * bulkmodulus[itype][itype]) / + (3.141592653589793 * cutsq[itype][jtype] * cutsq[itype][jtype]); + rk = (kshort * vfrac[j]) * (dr / cut[itype][jtype]); + + if (r > 0.0) fpair = -(rk/r); + else fpair = 0.0; + + f[i][0] += delx*fpair; + f[i][1] += dely*fpair; + f[i][2] += delz*fpair; + if (newton_pair || j < nlocal) { + f[j][0] -= delx*fpair; + f[j][1] -= dely*fpair; + f[j][2] -= delz*fpair; + } + + if (eflag) evdwl = 0.5*rk*dr; + if (evflag) ev_tally(i,j,nlocal,newton_pair,evdwl,0.0, + fpair*vfrac[i],delx,dely,delz); + } + } + } + + // grow bond forces array if necessary + + int maxpartner = 0; + for (i = 0; i < nlocal; i++) maxpartner = MAX(maxpartner,npartner[i]); + + + if (atom->nmax > nmax) { + memory->destroy(s0_new); + memory->destroy(theta); + nmax = atom->nmax; + memory->create(s0_new,nmax,"pair:s0_new"); + memory->create(theta,nmax,"pair:theta"); + + } + + // ******** temp array to store Plastic extension *********** /// + double deviatorPlasticExtTemp[nlocal][maxpartner]; + for (int ii = 0; ii < nlocal; ii++) { + for (int kk = 0; kk < maxpartner; kk++) { + deviatorPlasticExtTemp[ii][kk] = 0.0; + } + } + // ******** temp array to store Plastic extension *********** /// + + + + // compute the dilatation on each particle + + compute_dilatation(); + + // communicate dilatation (theta) of each particle + comm->forward_comm_pair(this); + + // communicate weighted volume (wvolume) upon every reneighbor + + if (neighbor->ago == 0) + comm->forward_comm_fix(modify->fix[ifix_peri]); + + // volume-dependent part of the energy + + if (eflag) { + for (i = 0; i < nlocal; i++) { + itype = type[i]; + if (eflag_global) + eng_vdwl += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]); + if (eflag_atom) + eatom[i] += 0.5 * bulkmodulus[itype][itype] * (theta[i] * theta[i]); + } + } + + // loop over my particles and their partners + // partner list contains all bond partners, so I-J appears twice + // if bond already broken, skip this partner + // first = true if this is first neighbor of particle i + + bool first; + double omega_minus, omega_plus, omega; + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + itype = type[i]; + jnum = npartner[i]; + first = true; + + + double yieldStress = m_yieldstress[itype][itype]; + double horizon = cut[itype][itype]; + double tdnorm = compute_DeviatoricForceStateNorm(i); + double pointwiseYieldvalue = 25.0 * yieldStress * + yieldStress / 8 / M_PI / pow(horizon,5); + + + double fsurf = (tdnorm * tdnorm)/2 - pointwiseYieldvalue; + bool elastic = true; + + double alphavalue = (15 * shearmodulus[itype][itype]) /wvolume[i]; + + + if (fsurf>0) { + elastic = false; + deltalambda = ((tdnorm /sqrt(2.0 * pointwiseYieldvalue)) - 1.0) / alphavalue; + double templambda = lambdaValue[i]; + lambdaValue[i] = templambda + deltalambda; + } + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + j = atom->map(partner[i][jj]); + // check if lost a partner without first breaking bond + + if (j < 0) { + partner[i][jj] = 0; + continue; + } + + // compute force density, add to PD equation of motion + + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (periodic) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + jtype = type[j]; + delta = cut[itype][jtype]; + r = sqrt(rsq); + dr = r - r0[i][jj]; + + // avoid roundoff errors + + if (fabs(dr) < 2.2204e-016) { + dr = 0.0; + } + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + omega_minus = influence_function(delx0,dely0,delz0); + + //Elastic Part + rk = ((3.0 * bulkmodulus[itype][itype]) * ( (omega_plus * theta[i] / wvolume[i]) + + ( omega_minus * theta[j] / wvolume[j] ) ) ) * r0[i][jj]; + + if (r > 0.0) fbond = -((rk/r) * vfrac[j] * vfrac_scale); + else fbond = 0.0; + + //Plastic part + + double deviatoric_extension = dr - (theta[i]* r0[i][jj] / 3.0); + edpNp1 = deviatorPlasticextension[i][jj]; + + + double tdtrialValue = ( 15 * shearmodulus[itype][itype]) * + ( (omega_plus / wvolume[i]) + (omega_minus / wvolume[j]) ) * + (deviatoric_extension - edpNp1); + + if(elastic) { + rkNew = tdtrialValue; + } + else { + rkNew = (sqrt(2.0*pointwiseYieldvalue) * tdtrialValue) / tdnorm; + deviatorPlasticExtTemp[i][jj] = edpNp1 + rkNew * deltalambda; + } + + + if (r > 0.0) fbondElastoPlastic = -((rkNew/r) * vfrac[j] * vfrac_scale); + else fbondElastoPlastic = 0.0; + + + // total Force state: elastic + plastic + fbondFinal=fbond+fbondElastoPlastic; + fbond=fbondFinal; + + + f[i][0] += delx*fbond; + f[i][1] += dely*fbond; + f[i][2] += delz*fbond; + + + // since I-J is double counted, set newton off & use 1/2 factor and I,I + + if (eflag) evdwl = (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * + omega_plus * deviatoric_extension * deviatoric_extension) + + (0.5 * 15 * (shearmodulus[itype][itype]/wvolume[i]) * + omega_plus * (deviatoric_extension - edpNp1) * + (deviatoric_extension-edpNp1)) * vfrac[j] * vfrac_scale; + if (evflag) ev_tally(i,i,nlocal,0,0.5*evdwl,0.0, + 0.5*fbond*vfrac[i],delx,dely,delz); + + // find stretch in bond I-J and break if necessary + // use s0 from previous timestep + + + stretch = dr / r0[i][jj]; + if (stretch > MIN(s0[i],s0[j])) partner[i][jj] = 0; + + // update s0 for next timestep + + if (first) + s0_new[i] = s00[itype][jtype] - (alpha[itype][jtype] * stretch); + else + s0_new[i] = MAX(s0_new[i],s00[itype][jtype] - + (alpha[itype][jtype] * stretch)); + + first = false; + } + } + + // store new s0 + + for (i = 0; i < nlocal; i++) s0[i] = s0_new[i]; + + for (i = 0; i < nlocal; i++) { + jnum = npartner[i]; + for (jj = 0; jj < jnum; jj++) { + double temp_data = deviatorPlasticExtTemp[i][jj]; + deviatorPlasticextension[i][jj] = temp_data; + } + } + +} + +/* ---------------------------------------------------------------------- + allocate all arrays +------------------------------------------------------------------------- */ + +void PairPeriEPS::allocate() +{ + allocated = 1; + int n = atom->ntypes; + + memory->create(setflag,n+1,n+1,"pair:setflag"); + for (int i = 1; i <= n; i++) + for (int j = i; j <= n; j++) + setflag[i][j] = 0; + + memory->create(cutsq,n+1,n+1,"pair:cutsq"); + memory->create(bulkmodulus,n+1,n+1,"pair:bulkmodulus"); + memory->create(shearmodulus,n+1,n+1,"pair:shearmodulus"); + memory->create(s00,n+1,n+1,"pair:s00"); + memory->create(alpha,n+1,n+1,"pair:alpha"); + memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(m_yieldstress,n+1,n+1,"pair:m_yieldstress"); +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairPeriEPS::settings(int narg, char **arg) +{ + if (narg) error->all(FLERR,"Illegal pair_style command"); +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairPeriEPS::coeff(int narg, char **arg) +{ + if (narg != 8) error->all(FLERR,"Incorrect args for pair coefficients"); + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(arg[0],atom->ntypes,ilo,ihi); + force->bounds(arg[1],atom->ntypes,jlo,jhi); + + double bulkmodulus_one = atof(arg[2]); + double shearmodulus_one = atof(arg[3]); + double cut_one = atof(arg[4]); + double s00_one = atof(arg[5]); + double alpha_one = atof(arg[6]); + double myieldstress_one = atof(arg[7]); + + int count = 0; + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + bulkmodulus[i][j] = bulkmodulus_one; + shearmodulus[i][j] = shearmodulus_one; + cut[i][j] = cut_one; + s00[i][j] = s00_one; + alpha[i][j] = alpha_one; + m_yieldstress[i][j] = myieldstress_one; + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairPeriEPS::init_one(int i, int j) +{ + if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set"); + + bulkmodulus[j][i] = bulkmodulus[i][j]; + shearmodulus[j][i] = shearmodulus[i][j]; + s00[j][i] = s00[i][j]; + alpha[j][i] = alpha[i][j]; + cut[j][i] = cut[i][j]; + m_yieldstress[j][i] = m_yieldstress[i][j]; + return cut[i][j]; +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairPeriEPS::init_style() +{ + // error checks + + if (!atom->peri_flag) + error->all(FLERR,"Pair style peri requires atom style peri"); + if (atom->map_style == 0) + error->all(FLERR,"Pair peri requires an atom map, see atom_modify"); + + if (domain->lattice == NULL) + error->all(FLERR,"Pair peri requires a lattice be defined"); + if (domain->lattice->xlattice != domain->lattice->ylattice || + domain->lattice->xlattice != domain->lattice->zlattice || + domain->lattice->ylattice != domain->lattice->zlattice) + error->all(FLERR,"Pair peri lattice is not identical in x, y, and z"); + + // if first init, create Fix needed for storing fixed neighbors + + if (ifix_peri == -1) { + char **fixarg = new char*[3]; + fixarg[0] = (char *) "PERI_NEIGH"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "PERI_NEIGH"; + modify->add_fix(3,fixarg); + delete [] fixarg; + } + + // find associated PERI_NEIGH fix that must exist + // could have changed locations in fix list since created + + for (int i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"PERI_NEIGH") == 0) ifix_peri = i; + if (ifix_peri == -1) error->all(FLERR,"Fix peri neigh does not exist"); + + neighbor->request(this); +} + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file +------------------------------------------------------------------------- */ + +void PairPeriEPS::write_restart(FILE *fp) +{ + int i,j; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + fwrite(&setflag[i][j],sizeof(int),1,fp); + if (setflag[i][j]) { + fwrite(&bulkmodulus[i][j],sizeof(double),1,fp); + fwrite(&shearmodulus[i][j],sizeof(double),1,fp); + fwrite(&s00[i][j],sizeof(double),1,fp); + fwrite(&alpha[i][j],sizeof(double),1,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + fwrite(&m_yieldstress[i][j],sizeof(double),1,fp); + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts +------------------------------------------------------------------------- */ + +void PairPeriEPS::read_restart(FILE *fp) +{ + allocate(); + + int i,j; + int me = comm->me; + for (i = 1; i <= atom->ntypes; i++) + for (j = i; j <= atom->ntypes; j++) { + if (me == 0) fread(&setflag[i][j],sizeof(int),1,fp); + MPI_Bcast(&setflag[i][j],1,MPI_INT,0,world); + if (setflag[i][j]) { + if (me == 0) { + fread(&bulkmodulus[i][j],sizeof(double),1,fp); + fread(&shearmodulus[i][j],sizeof(double),1,fp); + fread(&s00[i][j],sizeof(double),1,fp); + fread(&alpha[i][j],sizeof(double),1,fp); + fread(&cut[i][j],sizeof(double),1,fp); + fread(&m_yieldstress[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&bulkmodulus[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&shearmodulus[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&s00[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&alpha[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + MPI_Bcast(&m_yieldstress[i][j],1,MPI_DOUBLE,0,world); + } + } +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays +------------------------------------------------------------------------- */ + +double PairPeriEPS::memory_usage() +{ + double bytes = 2 * nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + influence function definition +------------------------------------------------------------------------- */ + +double PairPeriEPS::influence_function(double xi_x, double xi_y, double xi_z) +{ + double r = sqrt(xi_x*xi_x + xi_y*xi_y + xi_z*xi_z); + double omega; + + if (fabs(r) < 2.2204e-016) + error->one(FLERR,"Divide by 0 in influence function"); + omega = 1.0/r; + return omega; +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriEPS::compute_dilatation() +{ + int i,j,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0; + double rsq,r,dr; + double delta; + + double **x = atom->x; + int *type = atom->type; + double **x0 = atom->x0; + int nlocal = atom->nlocal; + double *vfrac = atom->vfrac; + double vfrac_scale = 1.0; + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + + int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; + + // compute the dilatation theta + + for (i = 0; i < nlocal; i++) { + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + jnum = npartner[i]; + theta[i] = 0.0; + itype = type[i]; + + for (jj = 0; jj < jnum; jj++) { + + // if bond already broken, skip this partner + if (partner[i][jj] == 0) continue; + + // look up local index of this partner particle + j = atom->map(partner[i][jj]); + + // skip if particle is "lost" + if (j < 0) continue; + + // compute force density and add to PD equation of motion + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (periodic) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + + r = sqrt(rsq); + dr = r - r0[i][jj]; + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + jtype = type[j]; + delta = cut[itype][jtype]; + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + theta[i] += influence_function(delx0, dely0, delz0) * r0[i][jj] * dr * + vfrac[j] * vfrac_scale; + + } + + // if wvolume[i] is zero, then particle i has no bonds + // therefore, the dilatation is set to + + if (wvolume[i] != 0.0) theta[i] = (3.0/wvolume[i]) * theta[i]; + else theta[i] = 0; + } +} + +/* ---------------------------------------------------------------------- */ + +double PairPeriEPS::compute_DeviatoricForceStateNorm(int i) +{ + int j,jj,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz; + double xtmp0,ytmp0,ztmp0,delx0,dely0,delz0; + double rsq,r,dr; + double delta; + double tdtrial; + double norm = 0.0; + + double **x = atom->x; + int *type = atom->type; + double **x0 = atom->x0; + double *s0 = atom->s0; + int nlocal = atom->nlocal; + double *vfrac = atom->vfrac; + double vfrac_scale = 1.0; + + double lc = domain->lattice->xlattice; + double half_lc = 0.5*lc; + + double **r0 = ((FixPeriNeigh *) modify->fix[ifix_peri])->r0; + tagint **partner = ((FixPeriNeigh *) modify->fix[ifix_peri])->partner; + int *npartner = ((FixPeriNeigh *) modify->fix[ifix_peri])->npartner; + double *wvolume = ((FixPeriNeigh *) modify->fix[ifix_peri])->wvolume; + double **deviatorPlasticextension = + ((FixPeriNeigh *) modify->fix[ifix_peri])->deviatorPlasticextension; + + int periodic = domain->xperiodic || domain->yperiodic || domain->zperiodic; + + // compute the dilatation theta + + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + xtmp0 = x0[i][0]; + ytmp0 = x0[i][1]; + ztmp0 = x0[i][2]; + jnum = npartner[i]; + itype = type[i]; + + for (jj = 0; jj < jnum; jj++) { + if (partner[i][jj] == 0) continue; + j = atom->map(partner[i][jj]); + // check if lost a partner without first breaking bond + if (j < 0) { + partner[i][jj] = 0; + continue; + } + delx = xtmp - x[j][0]; + dely = ytmp - x[j][1]; + delz = ztmp - x[j][2]; + if (periodic) domain->minimum_image(delx,dely,delz); + rsq = delx*delx + dely*dely + delz*delz; + delx0 = xtmp0 - x0[j][0]; + dely0 = ytmp0 - x0[j][1]; + delz0 = ztmp0 - x0[j][2]; + if (periodic) domain->minimum_image(delx0,dely0,delz0); + r = sqrt(rsq); + dr = r - r0[i][jj]; + if (fabs(dr) < 2.2204e-016) dr = 0.0; + + // scale vfrac[j] if particle j near the horizon + double vfrac_scale; + + jtype = type[j]; + double delta = cut[itype][jtype]; + + // scale vfrac[j] if particle j near the horizon + + if ((fabs(r0[i][jj] - delta)) <= half_lc) + vfrac_scale = (-1.0/(2*half_lc))*(r0[i][jj]) + + (1.0 + ((delta - half_lc)/(2*half_lc) ) ); + else vfrac_scale = 1.0; + + double ed = dr - (theta[i] * r0[i][jj])/3; + double edPNP1 = deviatorPlasticextension[i][jj]; + + jtype = type[j]; + delta = cut[itype][jtype]; + + double omega_plus = influence_function(-1.0*delx0,-1.0*dely0,-1.0*delz0); + double omega_minus = influence_function(delx0,dely0,delz0); + + double stretch = dr / r0[i][jj]; + + tdtrial = ( 15 * shearmodulus[itype][itype]) * + ((omega_plus * theta[i] / wvolume[i]) + + ( omega_minus * theta[j] / wvolume[j] ) ) * (ed - edPNP1); + + norm += tdtrial * tdtrial * vfrac[j] * vfrac_scale; + } + return sqrt(norm); +} + + +/* ---------------------------------------------------------------------- + communication routines +---------------------------------------------------------------------- */ + +int PairPeriEPS::pack_comm(int n, int *list, double *buf, + int pbc_flag, int *pbc) +{ + int i,j,m; + + m = 0; + for (i = 0; i < n; i++) { + j = list[i]; + buf[m++] = theta[j]; + } + return 1; +} + +/* ---------------------------------------------------------------------- */ + +void PairPeriEPS::unpack_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) { + theta[i] = buf[m++]; + } +} diff --git a/src/PERI/pair_peri_eps.h b/src/PERI/pair_peri_eps.h new file mode 100644 index 0000000000..c6eacce1c6 --- /dev/null +++ b/src/PERI/pair_peri_eps.h @@ -0,0 +1,111 @@ +/* ---------------------------------------------------------------------- + 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. +------------------------------------------------------------------------- */ + +#ifdef PAIR_CLASS + +PairStyle(peri/eps,PairPeriEPS) + +#else + +#ifndef LMP_PAIR_PERI_EPS_H +#define LMP_PAIR_PERI_EPS_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairPeriEPS : public Pair { + public: + double *theta; + + PairPeriEPS(class LAMMPS *); + virtual ~PairPeriEPS(); + int pack_comm(int, int *, double *, int, int *); + void unpack_comm(int, int, double *); + + virtual void compute(int, int); + void settings(int, char **); + void coeff(int, char **); + double init_one(int, int); + void init_style(); + void write_restart(FILE *); + void read_restart(FILE *); + void write_restart_settings(FILE *) {} + void read_restart_settings(FILE *) {} + double memory_usage(); + double influence_function(double, double, double); + void compute_dilatation(); + double compute_DeviatoricForceStateNorm(int); + + protected: + int ifix_peri; + double **bulkmodulus; + double **shearmodulus; + double **s00, **alpha; + double **cut, **m_yieldstress; //NEW: **m_yieldstress + + double *s0_new; + int nmax; + + void allocate(); +}; + +} + +#endif +#endif + +/* ERROR/WARNING messages: + +E: Illegal ... command + +Self-explanatory. Check the input script syntax and compare to the +documentation for the command. You can use -echo screen as a +command-line option when running LAMMPS to see the offending line. + +E: Incorrect args for pair coefficients + +Self-explanatory. Check the input script or data file. + +E: All pair coeffs are not set + +All pair coefficients must be set in the data file or by the +pair_coeff command before running a simulation. + +E: Pair style peri requires atom style peri + +Self-explanatory. + +E: Pair peri requires an atom map, see atom_modify + +Even for atomic systems, an atom map is required to find Peridynamic +bonds. Use the atom_modify command to define one. + +E: Pair peri requires a lattice be defined + +Use the lattice command for this purpose. + +E: Pair peri lattice is not identical in x, y, and z + +The lattice defined by the lattice command must be cubic. + +E: Fix peri neigh does not exist + +Somehow a fix that the pair style defines has been deleted. + +E: Divide by 0 in influence function of pair peri/lps + +This should not normally occur. It is likely a problem with your +model. + +*/ diff --git a/src/PERI/pair_peri_lps.cpp b/src/PERI/pair_peri_lps.cpp index 639ac6fb43..c8ee9548a5 100644 --- a/src/PERI/pair_peri_lps.cpp +++ b/src/PERI/pair_peri_lps.cpp @@ -620,12 +620,11 @@ void PairPeriLPS::compute_dilatation() /* ---------------------------------------------------------------------- communication routines - ---------------------------------------------------------------------- */ + ---------------------------------------------------------------------- */ int PairPeriLPS::pack_comm(int n, int *list, double *buf, int pbc_flag, int *pbc) { - int i,j,m; m = 0; diff --git a/src/PERI/pair_peri_lps.h b/src/PERI/pair_peri_lps.h index c9b063feea..c3593e57ac 100644 --- a/src/PERI/pair_peri_lps.h +++ b/src/PERI/pair_peri_lps.h @@ -26,6 +26,8 @@ namespace LAMMPS_NS { class PairPeriLPS : public Pair { public: + double *theta; + PairPeriLPS(class LAMMPS *); virtual ~PairPeriLPS(); int pack_comm(int, int *, double *, int, int *); @@ -48,11 +50,10 @@ class PairPeriLPS : public Pair { int ifix_peri; double **bulkmodulus; double **shearmodulus; - double **s00, **alpha; + double **s00,**alpha; double **cut; double *s0_new; - double *theta; int nmax; void allocate(); diff --git a/src/PERI/pair_peri_ves.cpp b/src/PERI/pair_peri_ves.cpp index ec7852feb5..baa826ea82 100644 --- a/src/PERI/pair_peri_ves.cpp +++ b/src/PERI/pair_peri_ves.cpp @@ -260,7 +260,7 @@ void PairPeriVES::compute(int eflag, int vflag) itype = type[i]; jnum = npartner[i]; first = true; - + for (jj = 0; jj < jnum; jj++) { if (partner[i][jj] == 0) continue; j = atom->map(partner[i][jj]); diff --git a/src/PERI/pair_peri_ves.h b/src/PERI/pair_peri_ves.h index c21e60c9b8..fbf6ab8bd9 100644 --- a/src/PERI/pair_peri_ves.h +++ b/src/PERI/pair_peri_ves.h @@ -26,6 +26,8 @@ namespace LAMMPS_NS { class PairPeriVES : public Pair { public: + double *theta; + PairPeriVES(class LAMMPS *); virtual ~PairPeriVES(); int pack_comm(int, int *, double *, int, int *); @@ -49,10 +51,11 @@ class PairPeriVES : public Pair { double **bulkmodulus; double **shearmodulus; double **s00, **alpha; - double **cut, **m_lambdai, **m_taubi; //NEW: **m_lambdai,**m_taubi + double **cut; + double **m_lambdai; + double **m_taubi; double *s0_new; - double *theta; int nmax; void allocate();