From 29dcdec8756bf714fc76b6e35026082494bea4d3 Mon Sep 17 00:00:00 2001 From: Dan Stefan Bolintineanu Date: Thu, 10 Jan 2019 16:53:50 -0700 Subject: [PATCH] Separated templated pair granular from pair granular/multi --- src/GRANULAR/pair_granular.cpp | 1069 ++++++++++++--- src/GRANULAR/pair_granular.h | 33 +- src/GRANULAR/pair_granular_multi.cpp | 1837 ++++++++++++++++++++++++++ src/GRANULAR/pair_granular_multi.h | 108 ++ src/fix_neigh_history.cpp | 8 +- src/pair.h | 4 + 6 files changed, 2850 insertions(+), 209 deletions(-) create mode 100644 src/GRANULAR/pair_granular_multi.cpp create mode 100644 src/GRANULAR/pair_granular_multi.h diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 9ee78686c6..c75c80fea5 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -48,9 +48,9 @@ using namespace MathConst; #define EPSILON 1e-10 -enum {VELOCITY, VISCOELASTIC, TSUJI}; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; -enum {TANGENTIAL_MINDLIN, TANGENTIAL_NOHISTORY}; +enum {VELOCITY, VISCOELASTIC, TSUJI}; +enum {TANGENTIAL_NOHISTORY, TANGENTIAL_MINDLIN}; enum {TWIST_NONE, TWIST_NOHISTORY, TWIST_SDS, TWIST_MARSHALL}; enum {ROLL_NONE, ROLL_NOHISTORY, ROLL_SDS}; @@ -105,12 +105,6 @@ PairGranular::~PairGranular() memory->destroy(roll_coeffs); memory->destroy(twist_coeffs); - memory->destroy(normal); - memory->destroy(damping); - memory->destroy(tangential); - memory->destroy(roll); - memory->destroy(twist); - delete [] onerad_dynamic; delete [] onerad_frozen; delete [] maxrad_dynamic; @@ -119,12 +113,730 @@ PairGranular::~PairGranular() memory->destroy(mass_rigid); } -void PairGranular::compute(int eflag, int vflag) +void PairGranular::compute(int eflag, int vflag){ +#ifdef TEMPLATED_PAIR_GRANULAR + if (normal == HOOKE){ + if (damping == VELOCITY){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<0,0,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<0,0,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<0,0,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<0,0,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<0,0,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<0,0,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<0,0,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<0,0,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,0,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,0,1,3,2>(eflag, vflag); + } + } + } + else if (damping == VISCOELASTIC){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<0,1,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<0,1,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<0,1,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<0,1,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<0,1,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<0,1,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<0,1,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<0,1,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,1,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,1,1,3,2>(eflag, vflag); + } + } + } + else if (damping == TSUJI){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<0,2,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<0,2,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<0,2,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<0,2,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<0,2,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<0,2,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<0,2,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<0,2,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<0,2,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<0,2,1,3,2>(eflag, vflag); + } + } + } + } + else if (normal == HERTZ){ + if (damping == VELOCITY){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<1,0,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<1,0,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<1,0,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<1,0,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<1,0,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<1,0,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<1,0,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<1,0,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,0,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,0,1,3,2>(eflag, vflag); + } + } + } + else if (damping == VISCOELASTIC){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<1,1,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<1,1,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<1,1,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<1,1,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<1,1,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<1,1,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<1,1,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<1,1,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,1,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,1,1,3,2>(eflag, vflag); + } + } + } + else if (damping == TSUJI){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<1,2,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<1,2,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<1,2,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<1,2,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<1,2,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<1,2,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<1,2,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<1,2,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<1,2,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<1,2,1,3,2>(eflag, vflag); + } + } + } + } + else if (normal == HERTZ_MATERIAL){ + if (damping == VELOCITY){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<2,0,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<2,0,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<2,0,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<2,0,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<2,0,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<2,0,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<2,0,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<2,0,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,0,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,0,1,3,2>(eflag, vflag); + } + } + } + else if (damping == VISCOELASTIC){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<2,1,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<2,1,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<2,1,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<2,1,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<2,1,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<2,1,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<2,1,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<2,1,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,1,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,1,1,3,2>(eflag, vflag); + } + } + } + else if (damping == TSUJI){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<2,2,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<2,2,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<2,2,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<2,2,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<2,2,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<2,2,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<2,2,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<2,2,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<2,2,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<2,2,1,3,2>(eflag, vflag); + } + } + } + } + else if (normal == DMT){ + if (damping == VELOCITY){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<3,0,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<3,0,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<3,0,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<3,0,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<3,0,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<3,0,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<3,0,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<3,0,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,0,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,0,1,3,2>(eflag, vflag); + } + } + } + else if (damping == VISCOELASTIC){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<3,1,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<3,1,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<3,1,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<3,1,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<3,1,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<3,1,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<3,1,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<3,1,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,1,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,1,1,3,2>(eflag, vflag); + } + } + } + else if (damping == TSUJI){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<3,2,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<3,2,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<3,2,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<3,2,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<3,2,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<3,2,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<3,2,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<3,2,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<3,2,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<3,2,1,3,2>(eflag, vflag); + } + } + } + } + else if (normal == JKR){ + if (damping == VELOCITY){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<4,0,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<4,0,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<4,0,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<4,0,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<4,0,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<4,0,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<4,0,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<4,0,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,0,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,0,1,3,2>(eflag, vflag); + } + } + } + else if (damping == VISCOELASTIC){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<4,1,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<4,1,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<4,1,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<4,1,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<4,1,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<4,1,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<4,1,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<4,1,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,1,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,1,1,3,2>(eflag, vflag); + } + } + } + else if (damping == TSUJI){ + if (tangential == TANGENTIAL_NOHISTORY){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<4,2,0,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,0,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,0,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<4,2,0,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,0,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,0,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<4,2,0,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,0,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,0,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<4,2,0,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,0,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,0,3,2>(eflag, vflag); + } + } + else if (tangential == TANGENTIAL_MINDLIN){ + if (twist == TWIST_NONE){ + if (roll == ROLL_NONE) compute_templated<4,2,1,0,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,1,0,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,1,0,2>(eflag, vflag); + } + else if (twist == TWIST_NOHISTORY){ + if (roll == ROLL_NONE) compute_templated<4,2,1,1,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,1,1,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,1,1,2>(eflag, vflag); + } + else if (twist == TWIST_SDS){ + if (roll == ROLL_NONE) compute_templated<4,2,1,2,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,1,2,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,1,2,2>(eflag, vflag); + } + else if (twist == TWIST_MARSHALL){ + if (roll == ROLL_NONE) compute_templated<4,2,1,3,0>(eflag, vflag); + else if (roll == ROLL_NOHISTORY) compute_templated<4,2,1,3,1>(eflag, vflag); + else if (roll == ROLL_SDS) compute_templated<4,2,1,3,2>(eflag, vflag); + } + } + } + } + +#else + compute_untemplated(Tp_normal, Tp_damping, Tp_tangential, + Tp_roll, Tp_twist, + eflag, vflag); +#endif +} + +#ifdef TEMPLATED_PAIR_GRANULAR +template < int Tp_normal, int Tp_damping, int Tp_tangential, + int Tp_roll, int Tp_twist > +void PairGranular::compute_templated(int eflag, int vflag) +#else +void PairGranular::compute_untemplated + (int Tp_normal, int Tp_damping, int Tp_tangential, + int Tp_roll, int Tp_twist, int eflag, int vflag) +#endif { int i,j,ii,jj,inum,jnum,itype,jtype; double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; double radi,radj,radsum,rsq,r,rinv,rsqinv; - double Reff, delta, dR, dR2, sqdR; + double Reff, delta, dR, dR2; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; double wr1,wr2,wr3; @@ -235,21 +947,21 @@ void PairGranular::compute(int eflag, int vflag) Reff = radi*radj/(radi+radj); touchflag = false; - if (normal[itype][jtype] == JKR){ + if (Tp_normal == JKR){ if (touch[jj]){ R2 = Reff*Reff; coh = normal_coeffs[itype][jtype][3]; a = cbrt(9.0*M_PI*coh*R2/(4*E)); delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E); dist_pulloff = radsum-delta_pulloff; - touchflag = (rsq <= dist_pulloff*dist_pulloff); + touchflag = (rsq < dist_pulloff*dist_pulloff); } else{ - touchflag = (rsq <= radsum*radsum); + touchflag = (rsq < radsum*radsum); } } else{ - touchflag = (rsq <= radsum*radsum); + touchflag = (rsq < radsum*radsum); } if (!touchflag){ @@ -296,8 +1008,10 @@ void PairGranular::compute(int eflag, int vflag) delta = radsum - r; dR = delta*Reff; - if (normal[itype][jtype] == JKR){ + if (Tp_normal == JKR){ touch[jj] = 1; + R2=Reff*Reff; + coh = normal_coeffs[itype][jtype][3]; dR2 = dR*dR; t0 = coh*coh*R2*R2*E; t1 = PI27SQ*t0; @@ -317,22 +1031,22 @@ void PairGranular::compute(int eflag, int vflag) else{ knfac = E; //Hooke Fne = knfac*delta; - if (normal[itype][jtype] != HOOKE) - a = sqdR = sqrt(dR); + if (Tp_normal != HOOKE) + a = sqrt(dR); Fne *= a; - if (normal[itype][jtype] == DMT) + if (Tp_normal == DMT) Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } //Consider restricting Hooke to only have 'velocity' as an option for damping? - if (damping[itype][jtype] == VELOCITY){ + if (Tp_damping == VELOCITY){ damp_normal = normal_coeffs[itype][jtype][1]; } - else if (damping[itype][jtype] == VISCOELASTIC){ - if (normal[itype][jtype] == HOOKE) a = sqdR = sqrt(dR); - damp_normal = normal_coeffs[itype][jtype][1]*sqdR*meff; + else if (Tp_damping == VISCOELASTIC){ + if (Tp_normal == HOOKE) a = sqrt(dR); + damp_normal = normal_coeffs[itype][jtype][1]*a*meff; } - else if (damping[itype][jtype] == TSUJI){ + else if (Tp_damping == TSUJI){ damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); } @@ -367,11 +1081,14 @@ void PairGranular::compute(int eflag, int vflag) history = &allhistory[size_history*jj]; } - Fcrit = fabs(Fne); - if (normal[itype][jtype] == JKR){ + + if (Tp_normal == JKR){ F_pulloff = 3*M_PI*coh*Reff; Fcrit = fabs(Fne + 2*F_pulloff); } + else{ + Fcrit = fabs(Fne); + } //------------------------------ //Tangential forces @@ -379,7 +1096,7 @@ void PairGranular::compute(int eflag, int vflag) k_tangential = tangential_coeffs[itype][jtype][0]; damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; - if (tangential_history){ + if (Tp_tangential > 0){ shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -436,7 +1153,7 @@ void PairGranular::compute(int eflag, int vflag) // Rolling resistance //**************************************** - if (roll[itype][jtype] != ROLL_NONE){ + if (Tp_roll != ROLL_NONE){ relrot1 = omega[i][0] - omega[j][0]; relrot2 = omega[i][1] - omega[j][1]; relrot3 = omega[i][2] - omega[j][2]; @@ -451,7 +1168,7 @@ void PairGranular::compute(int eflag, int vflag) if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; else vrlmaginv = 0.0; - if (roll_history){ + if (Tp_roll > 1){ int rhist0 = roll_history_index; int rhist1 = rhist0 + 1; int rhist2 = rhist1 + 1; @@ -515,9 +1232,9 @@ void PairGranular::compute(int eflag, int vflag) //**************************************** // Twisting torque, including history effects //**************************************** - if (twist[itype][jtype] != TWIST_NONE){ + if (Tp_twist != TWIST_NONE){ magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (twist[itype][jtype] == TWIST_MARSHALL){ + if (Tp_twist == TWIST_MARSHALL){ k_twist = 0.5*k_tangential*a*a;; //eq 32 damp_twist = 0.5*damp_tangential*a*a; mu_twist = TWOTHIRDS*a; @@ -527,7 +1244,7 @@ void PairGranular::compute(int eflag, int vflag) damp_twist = twist_coeffs[itype][jtype][1]; mu_twist = twist_coeffs[itype][jtype][2]; } - if (twist_history){ + if (Tp_twist > 1){ if (historyupdate){ history[twist_history_index] += magtwist*dt; } @@ -562,7 +1279,7 @@ void PairGranular::compute(int eflag, int vflag) torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; - if (twist[itype][jtype] != TWIST_NONE){ + if (Tp_twist != TWIST_NONE){ tortwist1 = magtortwist * nx; tortwist2 = magtortwist * ny; tortwist3 = magtortwist * nz; @@ -572,7 +1289,7 @@ void PairGranular::compute(int eflag, int vflag) torque[i][2] += tortwist3; } - if (roll[itype][jtype] != ROLL_NONE){ + if (Tp_roll != ROLL_NONE){ torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr torroll2 = Reff*(nz*fr1 - nx*fr3); torroll3 = Reff*(nx*fr2 - ny*fr1); @@ -591,12 +1308,12 @@ void PairGranular::compute(int eflag, int vflag) torque[j][1] -= radj*tor2; torque[j][2] -= radj*tor3; - if (twist[itype][jtype] != TWIST_NONE){ + if (Tp_twist != TWIST_NONE){ torque[j][0] -= tortwist1; torque[j][1] -= tortwist2; torque[j][2] -= tortwist3; } - if (roll[itype][jtype] != ROLL_NONE){ + if (Tp_roll != ROLL_NONE){ torque[j][0] -= torroll1; torque[j][1] -= torroll2; torque[j][2] -= torroll3; @@ -631,12 +1348,6 @@ void PairGranular::allocate() memory->create(roll_coeffs,n+1,n+1,3,"pair:roll_coeffs"); memory->create(twist_coeffs,n+1,n+1,3,"pair:twist_coeffs"); - memory->create(normal,n+1,n+1,"pair:normal"); - memory->create(damping,n+1,n+1,"pair:damping"); - memory->create(tangential,n+1,n+1,"pair:tangential"); - memory->create(roll,n+1,n+1,"pair:roll"); - memory->create(twist,n+1,n+1,"pair:twist"); - onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; maxrad_dynamic = new double[n+1]; @@ -654,12 +1365,14 @@ void PairGranular::settings(int narg, char **arg) int iarg = 0; //Some defaults - normal_global = HERTZ; - damping_global = VISCOELASTIC; - tangential_global = TANGENTIAL_MINDLIN; - roll_global = ROLL_NONE; - twist_global = TWIST_NONE; - tangential_history = roll_history = twist_history = 0; + normal = HERTZ; + damping = VISCOELASTIC; + tangential = TANGENTIAL_MINDLIN; + roll = ROLL_NONE; + twist = TWIST_NONE; + + tangential_history = 1; + roll_history = twist_history = 0; int normal_set, tangential_set; normal_set = tangential_set = 0; @@ -667,8 +1380,7 @@ void PairGranular::settings(int narg, char **arg) while (iarg < narg){ if (strcmp(arg[iarg], "hooke") == 0){ if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hooke option"); - normal_global = HOOKE; - normal_set = 1; + normal = HOOKE; memory->create(normal_coeffs_global, 2, "pair:normal_coeffs_global"); normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping @@ -677,7 +1389,7 @@ void PairGranular::settings(int narg, char **arg) else if (strcmp(arg[iarg], "hertz") == 0){ int num_coeffs = 2; if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal_global = HERTZ; + normal = HERTZ; normal_set = 1; memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn @@ -687,7 +1399,7 @@ void PairGranular::settings(int narg, char **arg) else if (strcmp(arg[iarg], "hertz/material") == 0){ int num_coeffs = 3; if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal_global = HERTZ_MATERIAL; + normal = HERTZ_MATERIAL; normal_set = 1; memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E (Young's modulus) @@ -697,7 +1409,7 @@ void PairGranular::settings(int narg, char **arg) } else if (strcmp(arg[iarg], "dmt") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); - normal_global = DMT; + normal = DMT; normal_set = 1; memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //4/3 E @@ -709,7 +1421,7 @@ void PairGranular::settings(int narg, char **arg) else if (strcmp(arg[iarg], "jkr") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for JKR option"); beyond_contact = 1; - normal_global = JKR; + normal = JKR; normal_set = 1; memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //E @@ -719,25 +1431,25 @@ void PairGranular::settings(int narg, char **arg) iarg += 5; } else if (strcmp(arg[iarg], "damp_velocity") == 0){ - damping_global = VELOCITY; + damping = VELOCITY; iarg += 1; } else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ - damping_global = VISCOELASTIC; + damping = VISCOELASTIC; iarg += 1; } else if (strcmp(arg[iarg], "damp_tsuji") == 0){ - damping_global = TSUJI; + damping = TSUJI; iarg += 1; } else if (strstr(arg[iarg], "tangential") != NULL){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for tangential model"); if (strstr(arg[iarg+1], "nohistory") != NULL){ - tangential_global = TANGENTIAL_NOHISTORY; + tangential = TANGENTIAL_NOHISTORY; tangential_set = 1; } else if (strstr(arg[iarg+1], "mindlin") != NULL){ - tangential_global = TANGENTIAL_MINDLIN; + tangential = TANGENTIAL_MINDLIN; tangential_history = 1; tangential_set = 1; } @@ -752,16 +1464,16 @@ void PairGranular::settings(int narg, char **arg) } else if (strstr(arg[iarg], "roll") != NULL){ if (strstr(arg[iarg+1], "none") != NULL){ - roll_global = ROLL_NONE; + roll = ROLL_NONE; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for rolling model"); if (strstr(arg[iarg+1], "nohistory") != NULL){ - roll_global = ROLL_NOHISTORY; + roll = ROLL_NOHISTORY; } else if (strstr(arg[iarg+1], "sds") != NULL){ - roll_global = ROLL_SDS; + roll = ROLL_SDS; roll_history = 1; } else{ @@ -776,11 +1488,11 @@ void PairGranular::settings(int narg, char **arg) } else if (strstr(arg[iarg], "twist") != NULL){ if (strstr(arg[iarg+1], "none") != NULL){ - twist_global = TWIST_NONE; + twist = TWIST_NONE; iarg += 2; } else if (strstr(arg[iarg+1], "marshall") != NULL){ - twist_global = TWIST_MARSHALL; + twist = TWIST_MARSHALL; twist_history = 1; iarg += 2; } @@ -788,10 +1500,10 @@ void PairGranular::settings(int narg, char **arg) if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for twist model"); memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); //To be filled later if (strstr(arg[iarg+1], "nohistory") != NULL){ - twist_global = TWIST_NOHISTORY; + twist = TWIST_NOHISTORY; } else if (strstr(arg[iarg+1], "sds") != NULL){ - twist_global = TWIST_SDS; + twist = TWIST_SDS; twist_history = 1; } else{ @@ -813,46 +1525,35 @@ void PairGranular::settings(int narg, char **arg) // The reason for the current setup is to remain true to existing pair gran/hooke etc. syntax, // where coeffs are set in the pair_style command, and a pair_coeff * * command is issued. - //Other option is to have two pair styles, e.g. pair granular and pair granular/multi, - // where granular/multi allows per-type coefficients, pair granular does not (this would also - // allow minor speed-up by templating pair granular) allocate(); double damp; - for (int i = 1; i <= atom->ntypes; i++){ - normal[i][i] = normal_global; - damping[i][i] = damping_global; - tangential[i][i] = tangential_global; - roll[i][i] = roll_global; - twist[i][i] = twist_global; + if (damping == TSUJI){ + double cor = normal_coeffs_global[1]; + damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ + 27.467*pow(cor,4)-18.022*pow(cor,5)+ + 4.8218*pow(cor,6); + } + else damp = normal_coeffs_global[1]; + for (int i = 1; i <= atom->ntypes; i++){ if (normal_set){ - if (damping_global == TSUJI){ - double cor = normal_coeffs_global[1]; - damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ - 27.467*pow(cor,4)-18.022*pow(cor,5)+ - 4.8218*pow(cor,6); - } - else damp = normal_coeffs_global[1]; normal_coeffs[i][i][0] = normal_coeffs_global[0]; normal_coeffs[i][i][1] = damp; - if (normal[i][i] != HOOKE && normal[i][i] != HERTZ){ + if (normal != HOOKE && normal != HERTZ){ normal_coeffs[i][i][2] = normal_coeffs_global[2]; } - if ((normal_global == JKR) || (normal_global == DMT)) + if ((normal == JKR) || (normal == DMT)) normal_coeffs[i][i][3] = normal_coeffs_global[3]; } if(tangential_set){ - tangential[i][i] = tangential_global; for (int k = 0; k < 3; k++) tangential_coeffs[i][i][k] = tangential_coeffs_global[k]; } - roll[i][i] = roll_global; - if (roll_global != ROLL_NONE) + if (roll != ROLL_NONE) for (int k = 0; k < 3; k++) roll_coeffs[i][i][k] = roll_coeffs_global[k]; - twist[i][i] = twist_global; - if (twist_global != TWIST_NONE && twist_global != TWIST_MARSHALL) + if (twist != TWIST_NONE && twist != TWIST_MARSHALL) for (int k = 0; k < 3; k++) twist_coeffs[i][i][k] = twist_coeffs_global[k]; @@ -866,7 +1567,7 @@ void PairGranular::settings(int narg, char **arg) void PairGranular::coeff(int narg, char **arg) { - int normal_local, damping_local, tangential_local, roll_local, twist_local; + int normal_set, damping_set, tangential_set, roll_set, twist_set; double *normal_coeffs_local; double *tangential_coeffs_local; double *roll_coeffs_local; @@ -886,78 +1587,83 @@ void PairGranular::coeff(int narg, char **arg) force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); - normal_local = tangential_local = roll_local = twist_local = -1; - damping_local = -1; + normal_set = damping_set = tangential_set = roll_set = twist_set = 0; int iarg = 2; while (iarg < narg){ if (strcmp(arg[iarg], "hooke") == 0){ if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hooke option"); - normal_local = HOOKE; + if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_set = 1; iarg += 3; } else if (strcmp(arg[iarg], "hertz") == 0){ int num_coeffs = 2; if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_local = HERTZ; + if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_set = 1; iarg += num_coeffs+1; } else if (strcmp(arg[iarg], "hertz/material") == 0){ int num_coeffs = 3; if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_local = HERTZ; + if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G + normal_set = 1; iarg += num_coeffs+1; } else if (strcmp(arg[iarg], "dmt") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); - normal_local = DMT; + if (normal != DMT) if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion + normal_set = 1; iarg += 5; } else if (strcmp(arg[iarg], "jkr") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for JKR option"); beyond_contact = 1; - normal_local = JKR; + if (normal != JKR) if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be consistent"); normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //E normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + normal_set = 1; iarg += 5; } else if (strcmp(arg[iarg], "damp_velocity") == 0){ - damping_local = VELOCITY; + if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent"); iarg += 1; } else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ - damping_local = VISCOELASTIC; + if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent"); iarg += 1; } else if (strcmp(arg[iarg], "damp_tsuji") == 0){ - damping_local = TSUJI; + if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be consistent"); iarg += 1; } else if (strstr(arg[iarg], "tangential") != NULL){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); if (strstr(arg[iarg+1], "nohistory") != NULL){ - tangential_local = TANGENTIAL_NOHISTORY; + if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be consistent"); } else if (strstr(arg[iarg+1], "mindlin") != NULL){ - tangential_local = TANGENTIAL_MINDLIN; + if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be consistent");; tangential_history = 1; } else{ error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); } + tangential_set = 1; tangential_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt tangential_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. @@ -966,21 +1672,22 @@ void PairGranular::coeff(int narg, char **arg) else if (strstr(arg[iarg], "rolling") != NULL){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); if (strstr(arg[iarg+1], "none") != NULL){ - roll_local = ROLL_NONE; + if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent"); iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); if (strstr(arg[iarg+1], "nohistory") != NULL){ - roll_local = ROLL_NOHISTORY; + if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent"); } else if (strstr(arg[iarg+1], "sds") != NULL){ - roll_local = ROLL_SDS; + if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be consistent"); roll_history = 1; } else{ error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized"); } + roll_set =1 ; roll_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt roll_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat roll_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. @@ -990,26 +1697,27 @@ void PairGranular::coeff(int narg, char **arg) else if (strstr(arg[iarg], "twist") != NULL){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); if (strstr(arg[iarg+1], "none") != NULL){ - twist_local = TWIST_NONE; + if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); iarg += 2; } else if (strstr(arg[iarg+1], "marshall") != NULL){ - twist_local = TWIST_MARSHALL; + if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); twist_history = 1; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); if (strstr(arg[iarg+1], "nohistory") != NULL){ - twist_local = TWIST_NOHISTORY; + if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); } else if (strstr(arg[iarg+1], "sds") != NULL){ - twist_local = TWIST_SDS; + if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be consistent"); twist_history = 1; } else{ error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized"); } + twist_set = 1; twist_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt twist_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat twist_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. @@ -1021,55 +1729,39 @@ void PairGranular::coeff(int narg, char **arg) int count = 0; double damp; - if (damping_local >= 0){ - if (normal_local == -1) - error->all(FLERR, "Illegal pair_coeff command, must specify normal model when setting damping model"); - if (damping_local == TSUJI){ - double cor; - cor = normal_coeffs_local[1]; - damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ - 27.467*pow(cor,4)-18.022*pow(cor,5)+ - 4.8218*pow(cor,6); - } - else damp = normal_coeffs_local[1]; + if (damping == TSUJI){ + double cor; + cor = normal_coeffs_local[1]; + damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ + 27.467*pow(cor,4)-18.022*pow(cor,5)+ + 4.8218*pow(cor,6); } + else damp = normal_coeffs_local[1]; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - if (normal_local >= 0){ - normal[i][j] = normal_local; + if (normal_set){ normal_coeffs[i][j][0] = normal_coeffs_local[0]; - if (damping_local == -1){ - damp = normal_coeffs_global[1]; - } normal_coeffs[i][j][1] = damp; - if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; - if ((normal_local == JKR) || (normal_local == DMT)) + if (normal != HERTZ && normal != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; + if ((normal == JKR) || (normal == DMT)) normal_coeffs[i][j][3] = normal_coeffs_local[3]; } - if (damping_local >= 0){ - damping[i][j] = damping_local; - } - if (tangential_local >= 0){ - tangential[i][j] = tangential_local; + if (tangential_set){ for (int k = 0; k < 3; k++) tangential_coeffs[i][j][k] = tangential_coeffs_local[k]; } - if (roll_local >= 0){ - roll[i][j] = roll_local; - if (roll_local != ROLL_NONE) + if (roll_set){ + if (roll != ROLL_NONE) for (int k = 0; k < 3; k++) roll_coeffs[i][j][k] = roll_coeffs_local[k]; } - if (twist_local >= 0){ - twist[i][j] = twist_local; - if (twist_local != TWIST_NONE && twist_local != TWIST_MARSHALL) + if (twist_set){ + if (twist != TWIST_NONE && twist != TWIST_MARSHALL) for (int k = 0; k < 3; k++) twist_coeffs[i][j][k] = twist_coeffs_local[k]; } - - if (normal_local >= 0 && tangential_local >= 0) setflag[i][j] = 1; - + setflag[i][j] = 1; count++; } } @@ -1101,9 +1793,7 @@ void PairGranular::init_style() use_history = tangential_history || roll_history || twist_history; //For JKR, will need fix/neigh/history to keep track of touch arrays - for (int i = 1; i <= atom->ntypes; i++) - for (int j = 1; j <= atom->ntypes; j++) - if (normal[i][j] == JKR) use_history = 1; + if (normal == JKR) use_history = 1; size_history = 3*tangential_history + 3*roll_history + twist_history; @@ -1181,13 +1871,13 @@ void PairGranular::init_style() if (ipour >= 0) { itype = i; double radmax = *((double *) modify->fix[ipour]->extract("radius",itype)); - if (normal[itype][itype] == JKR) radmax = radmax + 0.5*pulloff_distance(radmax, itype); + if (normal == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); onerad_dynamic[i] = radmax; } if (idep >= 0) { itype = i; double radmax = *((double *) modify->fix[idep]->extract("radius",itype)); - if (normal[itype][itype] == JKR) radmax = radmax + 0.5*pulloff_distance(radmax, itype); + if (normal == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); onerad_dynamic[i] = radmax; } } @@ -1199,8 +1889,8 @@ void PairGranular::init_style() for (i = 0; i < nlocal; i++){ double radius_cut = radius[i]; - if (normal[type[i]][type[i]] == JKR){ - radius_cut = radius[i] + 0.5*pulloff_distance(radius[i], type[i]); + if (normal == JKR){ + radius_cut = radius[i] - 0.5*pulloff_distance(radius[i], type[i]); } if (mask[i] & freeze_group_bit){ onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius_cut); @@ -1232,18 +1922,8 @@ double PairGranular::init_one(int i, int j) { double cutoff; if (setflag[i][j] == 0) { - if ((normal[i] != normal[j]) || - (damping[i] != damping[j]) || - (tangential[i] != tangential[j]) || - (roll[i] != roll[j]) || - (twist[i] != twist[j])){ - char str[512]; - sprintf(str,"Granular pair style functional forms are different, cannot mix coefficients for types %d and %d. \nThis combination must be set explicitly via pair_coeff command.",i,j); - error->one(FLERR,str); - } - - if (normal[i][j] != HOOKE && normal[i][j] != HERTZ){ + if (normal != HOOKE && normal != HERTZ){ normal_coeffs[i][j][0] = mix_stiffnessE(normal_coeffs[i][i][0], normal_coeffs[j][j][0], normal_coeffs[i][i][2], normal_coeffs[j][j][2]); normal_coeffs[i][j][2] = mix_stiffnessG(normal_coeffs[i][i][0], normal_coeffs[j][j][0], @@ -1254,19 +1934,19 @@ double PairGranular::init_one(int i, int j) } normal_coeffs[i][j][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); - if ((normal[i][i] == JKR) || (normal[i][i] == DMT)) + if ((normal == JKR) || (normal == DMT)) normal_coeffs[i][j][3] = mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]); for (int k = 0; k < 3; k++) tangential_coeffs[i][j][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); - if (roll[i][i] != ROLL_NONE){ + if (roll != ROLL_NONE){ for (int k = 0; k < 3; k++) roll_coeffs[i][j][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); } - if (twist[i][i] != TWIST_NONE && twist[i][i] != TWIST_MARSHALL){ + if (twist != TWIST_NONE && twist != TWIST_MARSHALL){ for (int k = 0; k < 3; k++) twist_coeffs[i][j][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]); } @@ -1305,15 +1985,15 @@ double PairGranular::init_one(int i, int j) void PairGranular::write_restart(FILE *fp) { int i,j; + fwrite(&normal,sizeof(int),1,fp); + fwrite(&damping,sizeof(int),1,fp); + fwrite(&tangential,sizeof(int),1,fp); + fwrite(&roll,sizeof(int),1,fp); + fwrite(&twist,sizeof(int),1,fp); 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(&normal[i][j],sizeof(int),1,fp); - fwrite(&damping[i][j],sizeof(int),1,fp); - fwrite(&tangential[i][j],sizeof(int),1,fp); - fwrite(&roll[i][j],sizeof(int),1,fp); - fwrite(&twist[i][j],sizeof(int),1,fp); fwrite(&normal_coeffs[i][j],sizeof(double),4,fp); fwrite(&tangential_coeffs[i][j],sizeof(double),3,fp); fwrite(&roll_coeffs[i][j],sizeof(double),3,fp); @@ -1333,28 +2013,30 @@ void PairGranular::read_restart(FILE *fp) allocate(); int i,j; int me = comm->me; + if (me == 0){ + fread(&normal,sizeof(int),1,fp); + fread(&damping,sizeof(int),1,fp); + fread(&tangential,sizeof(int),1,fp); + fread(&roll,sizeof(int),1,fp); + fread(&twist,sizeof(int),1,fp); + } + MPI_Bcast(&normal,1,MPI_INT,0,world); + MPI_Bcast(&damping,1,MPI_INT,0,world); + MPI_Bcast(&tangential,1,MPI_INT,0,world); + MPI_Bcast(&roll,1,MPI_INT,0,world); + MPI_Bcast(&twist,1,MPI_INT,0,world); 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(&normal[i][j],sizeof(int),1,fp); - fread(&damping[i][j],sizeof(int),1,fp); - fread(&tangential[i][j],sizeof(int),1,fp); - fread(&roll[i][j],sizeof(int),1,fp); - fread(&twist[i][j],sizeof(int),1,fp); fread(&normal_coeffs[i][j],sizeof(double),4,fp); fread(&tangential_coeffs[i][j],sizeof(double),3,fp); fread(&roll_coeffs[i][j],sizeof(double),3,fp); fread(&twist_coeffs[i][j],sizeof(double),3,fp); fread(&cut[i][j],sizeof(double),1,fp); } - MPI_Bcast(&normal[i][j],1,MPI_INT,0,world); - MPI_Bcast(&damping[i][j],1,MPI_INT,0,world); - MPI_Bcast(&tangential[i][j],1,MPI_INT,0,world); - MPI_Bcast(&roll[i][j],1,MPI_INT,0,world); - MPI_Bcast(&twist[i][j],1,MPI_INT,0,world); MPI_Bcast(&normal_coeffs[i][j],4,MPI_DOUBLE,0,world); MPI_Bcast(&tangential_coeffs[i][j],3,MPI_DOUBLE,0,world); MPI_Bcast(&roll_coeffs[i][j],3,MPI_DOUBLE,0,world); @@ -1380,7 +2062,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, { double radi,radj,radsum; double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, Reff; - double dR, dR2, sqdR; + double dR, dR2; double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; double vtr1,vtr2,vtr3,vrel; double mi,mj,meff,damp,ccel,tor1,tor2,tor3; @@ -1421,7 +2103,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Reff = radi*radj/(radi+radj); bool touchflag; - if (normal[itype][jtype] == JKR){ + if (normal == JKR){ R2 = Reff*Reff; coh = normal_coeffs[itype][jtype][3]; a = cbrt(9.0*M_PI*coh*R2/(4*E)); @@ -1513,7 +2195,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, delta = radsum - r; dR = delta*Reff; - if (normal[itype][jtype] == JKR){ + if (normal == JKR){ dR2 = dR*dR; t0 = coh*coh*R2*R2*E; t1 = PI27SQ*t0; @@ -1529,29 +2211,26 @@ double PairGranular::single(int i, int j, int itype, int jtype, a2 = a*a; knfac = FOURTHIRDS*E*a; Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); - if (damping[itype][jtype] == VISCOELASTIC) sqdR = sqrt(dR); } else{ knfac = E; Fne = knfac*delta; - if (normal[itype][jtype] != HOOKE) - a = sqdR = sqrt(dR); + if (normal != HOOKE) + a = sqrt(dR); Fne *= a; - if (normal[itype][jtype] == DMT) + if (normal == DMT) Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; } //Consider restricting Hooke to only have 'velocity' as an option for damping? - if (damping[itype][jtype] == VELOCITY){ + if (damping == VELOCITY){ damp_normal = normal_coeffs[itype][jtype][1]; } - else if (damping[itype][jtype] == VISCOELASTIC){ - if (normal[itype][jtype] == HOOKE) a = sqdR = sqrt(dR); - - - damp_normal = normal_coeffs[itype][jtype][1]*sqdR*meff; + else if (damping == VISCOELASTIC){ + if (normal == HOOKE) a = sqrt(dR); + damp_normal = normal_coeffs[itype][jtype][1]*a*meff; } - else if (damping[itype][jtype] == TSUJI){ + else if (damping == TSUJI){ damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); } @@ -1594,7 +2273,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, vrel = sqrt(vrel); Fcrit = fabs(Fne); - if (normal[itype][jtype] == JKR){ + if (normal == JKR){ F_pulloff = 3*M_PI*coh*Reff; Fcrit = fabs(Fne + 2*F_pulloff); } @@ -1603,7 +2282,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, //Tangential forces //------------------------------ k_tangential = tangential_coeffs[itype][jtype][0]; - if (normal[itype][jtype] != HOOKE){ + if (normal != HOOKE){ k_tangential *= a; } damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; @@ -1644,7 +2323,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, // Rolling resistance //**************************************** - if (roll[itype][jtype] != ROLL_NONE){ + if (roll != ROLL_NONE){ relrot1 = omega[i][0] - omega[j][0]; relrot2 = omega[i][1] - omega[j][1]; relrot3 = omega[i][2] - omega[j][2]; @@ -1705,9 +2384,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, //**************************************** // Twisting torque, including history effects //**************************************** - if (twist[itype][jtype] != TWIST_NONE){ + if (twist != TWIST_NONE){ magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (twist[itype][jtype] == TWIST_MARSHALL){ + if (twist == TWIST_MARSHALL){ k_twist = 0.5*k_tangential*a*a;; //eq 32 damp_twist = 0.5*damp_tangential*a*a; mu_twist = TWOTHIRDS*a; diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index 0e6f28c514..897316c907 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -28,7 +28,19 @@ class PairGranular : public Pair { public: PairGranular(class LAMMPS *); virtual ~PairGranular(); - virtual void compute(int, int); + + void compute(int, int); + // comment next line to turn off templating +#define TEMPLATED_PAIR_GRANULAR +#ifdef TEMPLATED_PAIR_GRANULAR + template < int Tp_normal, int Tp_damping, int Tp_tangential, + int Tp_roll, int Tp_twist> + void compute_templated(int, int); +#else + void compute_untemplated(int, int, int, int, int, + int, int); +#endif + virtual void settings(int, char **); virtual void coeff(int, char **); void init_style(); @@ -61,28 +73,27 @@ public: int nmax; // allocated size of mass_rigid virtual void allocate(); - int beyond_contact; - int nondefault_history_transfer; private: int size_history; - //Per-type models - int **normal, **damping, **tangential, **roll, **twist; - - int normal_global, damping_global; - int tangential_global, roll_global, twist_global; + //Models + int normal, damping, tangential, roll, twist; + //History flags int tangential_history, roll_history, twist_history; - int tangential_history_index; - int roll_history_index; - int twist_history_index; + //Indices of history entries + int tangential_history_index, roll_history_index, twist_history_index; + + //Coefficients declared in pair style command, used as default unless + // overwritten in pair coeff command double *normal_coeffs_global; double *tangential_coeffs_global; double *roll_coeffs_global; double *twist_coeffs_global; + //Per-type coefficients declared in pair coeff command double ***normal_coeffs; double ***tangential_coeffs; double ***roll_coeffs; diff --git a/src/GRANULAR/pair_granular_multi.cpp b/src/GRANULAR/pair_granular_multi.cpp new file mode 100644 index 0000000000..11381444a2 --- /dev/null +++ b/src/GRANULAR/pair_granular_multi.cpp @@ -0,0 +1,1837 @@ +/* ---------------------------------------------------------------------- +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: +Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL) +Leo Silbert (SNL), Gary Grest (SNL) +----------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_granular_multi.h" +#include "atom.h" +#include "atom_vec.h" +#include "domain.h" +#include "force.h" +#include "update.h" +#include "modify.h" +#include "fix.h" +#include "fix_neigh_history.h" +#include "comm.h" +#include "neighbor.h" +#include "neigh_list.h" +#include "neigh_request.h" +#include "memory.h" +#include "error.h" +#include "math_const.h" + +using namespace LAMMPS_NS; +using namespace MathConst; + +#define PI27SQ 266.47931882941264802866 // 27*PI**2 +#define THREEROOT3 5.19615242270663202362 // 3*sqrt(3) +#define SIXROOT6 14.69693845669906728801 // 6*sqrt(6) +#define INVROOT6 0.40824829046386307274 // 1/sqrt(6) +#define FOURTHIRDS 1.333333333333333 // 4/3 +#define TWOPI 6.28318530717959 // 2*PI + +#define EPSILON 1e-10 + +enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; +enum {VELOCITY, VISCOELASTIC, TSUJI}; +enum {TANGENTIAL_NOHISTORY, TANGENTIAL_MINDLIN}; +enum {TWIST_NONE, TWIST_NOHISTORY, TWIST_SDS, TWIST_MARSHALL}; +enum {ROLL_NONE, ROLL_NOHISTORY, ROLL_SDS}; + +/* ---------------------------------------------------------------------- */ + +PairGranularMulti::PairGranularMulti(LAMMPS *lmp) : Pair(lmp) +{ + single_enable = 1; + no_virial_fdotr_compute = 1; + fix_history = NULL; + + single_extra = 9; + svector = new double[single_extra]; + + neighprev = 0; + + nmax = 0; + mass_rigid = NULL; + + onerad_dynamic = NULL; + onerad_frozen = NULL; + maxrad_dynamic = NULL; + maxrad_frozen = NULL; + + dt = update->dt; + + // set comm size needed by this Pair if used with fix rigid + + comm_forward = 1; + + use_history = 0; + beyond_contact = 0; + nondefault_history_transfer = 0; + tangential_history_index = 0; + roll_history_index = twist_history_index = 0; + +} + +/* ---------------------------------------------------------------------- */ +PairGranularMulti::~PairGranularMulti() +{ + delete [] svector; + if (fix_history) modify->delete_fix("NEIGH_HISTORY"); + + if (allocated) { + memory->destroy(setflag); + memory->destroy(cutsq); + memory->destroy(cut); + + memory->destroy(normal_coeffs); + memory->destroy(tangential_coeffs); + memory->destroy(roll_coeffs); + memory->destroy(twist_coeffs); + + memory->destroy(normal); + memory->destroy(damping); + memory->destroy(tangential); + memory->destroy(roll); + memory->destroy(twist); + + delete [] onerad_dynamic; + delete [] onerad_frozen; + delete [] maxrad_dynamic; + delete [] maxrad_frozen; + } + memory->destroy(mass_rigid); +} + +void PairGranularMulti::compute(int eflag, int vflag) +{ + int i,j,ii,jj,inum,jnum,itype,jtype; + double xtmp,ytmp,ztmp,delx,dely,delz,fx,fy,fz,nx,ny,nz; + double radi,radj,radsum,rsq,r,rinv,rsqinv; + double Reff, delta, dR, dR2; + + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3; + double wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; + + double knfac, damp_normal; + double k_tangential, damp_tangential; + double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; + double fs, fs1, fs2, fs3; + + //For JKR + double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; + double t0, t1, t2, t3, t4, t5, t6; + double sqrt1, sqrt2, sqrt3, sqrt4; + + double mi,mj,meff,damp,ccel,tor1,tor2,tor3; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; + + //Rolling + double k_roll, damp_roll; + double roll1, roll2, roll3, torroll1, torroll2, torroll3; + double rollmag, rolldotn, scalefac; + double fr, fr1, fr2, fr3; + + //Twisting + double k_twist, damp_twist, mu_twist; + double signtwist, magtwist, magtortwist, Mtcrit; + double tortwist1, tortwist2, tortwist3; + + double shrmag,rsht; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *history,*allhistory,**firsthistory; + + bool touchflag; + + if (eflag || vflag) ev_setup(eflag,vflag); + else evflag = vflag_fdotr = 0; + + int historyupdate = 1; + if (update->setupflag) historyupdate = 0; + + // update rigid body info for owned & ghost atoms if using FixRigid masses + // body[i] = which body atom I is in, -1 if none + // mass_body = mass of each rigid body + + if (fix_rigid && neighbor->ago == 0){ + int tmp; + int *body = (int *) fix_rigid->extract("body",tmp); + double *mass_body = (double *) fix_rigid->extract("masstotal",tmp); + if (atom->nmax > nmax) { + memory->destroy(mass_rigid); + nmax = atom->nmax; + memory->create(mass_rigid,nmax,"pair:mass_rigid"); + } + int nlocal = atom->nlocal; + for (i = 0; i < nlocal; i++) + if (body[i] >= 0) mass_rigid[i] = mass_body[body[i]]; + else mass_rigid[i] = 0.0; + comm->forward_comm_pair(this); + } + + double **x = atom->x; + double **v = atom->v; + double **f = atom->f; + int *type = atom->type; + double **omega = atom->omega; + double **torque = atom->torque; + double *radius = atom->radius; + double *rmass = atom->rmass; + int *mask = atom->mask; + int nlocal = atom->nlocal; + int newton_pair = force->newton_pair; + + inum = list->inum; + ilist = list->ilist; + numneigh = list->numneigh; + firstneigh = list->firstneigh; + firsttouch = fix_history->firstflag; + firsthistory = fix_history->firstvalue; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + itype = type[i]; + xtmp = x[i][0]; + ytmp = x[i][1]; + ztmp = x[i][2]; + itype = type[i]; + radi = radius[i]; + touch = firsttouch[i]; + allhistory = firsthistory[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]; + jtype = type[j]; + rsq = delx*delx + dely*dely + delz*delz; + radj = radius[j]; + radsum = radi + radj; + + E = normal_coeffs[itype][jtype][0]; + Reff = radi*radj/(radi+radj); + touchflag = false; + + if (normal[itype][jtype] == JKR){ + if (touch[jj]){ + R2 = Reff*Reff; + coh = normal_coeffs[itype][jtype][3]; + a = cbrt(9.0*M_PI*coh*R2/(4*E)); + delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E); + dist_pulloff = radsum-delta_pulloff; + touchflag = (rsq < dist_pulloff*dist_pulloff); + } + else{ + touchflag = (rsq < radsum*radsum); + } + } + else{ + touchflag = (rsq < radsum*radsum); + } + + if (!touchflag){ + // unset non-touching neighbors + touch[jj] = 0; + history = &allhistory[size_history*jj]; + for (int k = 0; k < size_history; k++) history[k] = 0.0; + } + else{ + r = sqrt(rsq); + rinv = 1.0/r; + + nx = delx*rinv; + ny = dely*rinv; + nz = delz*rinv; + + // relative translational velocity + + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + vnnr = vr1*nx + vr2*ny + vr3*nz; //v_R . n + vn1 = nx*vnnr; + vn2 = ny*vnnr; + vn3 = nz*vnnr; + + // meff = effective mass of pair of particles + // if I or J part of rigid body, use body mass + // if I or J is frozen, meff is other particle + + mi = rmass[i]; + mj = rmass[j]; + if (fix_rigid) { + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; + } + + meff = mi*mj / (mi+mj); + if (mask[i] & freeze_group_bit) meff = mj; + if (mask[j] & freeze_group_bit) meff = mi; + + delta = radsum - r; + dR = delta*Reff; + if (normal[itype][jtype] == JKR){ + touch[jj] = 1; + R2=Reff*Reff; + coh = normal_coeffs[itype][jtype][3]; + dR2 = dR*dR; + t0 = coh*coh*R2*R2*E; + t1 = PI27SQ*t0; + t2 = 8*dR*dR2*E*E*E; + t3 = 4*dR2*E; + sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues + t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); + t5 = t3/t4 + t4/E; + sqrt2 = MAX(0, 2*dR + t5); + t6 = sqrt(sqrt2); + sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); + a = INVROOT6*(t6 + sqrt(sqrt3)); + a2 = a*a; + knfac = FOURTHIRDS*E*a; + Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); + } + else{ + knfac = E; //Hooke + Fne = knfac*delta; + if (normal[itype][jtype] != HOOKE) + a = sqrt(dR); + Fne *= a; + if (normal[itype][jtype] == DMT) + Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; + } + + //Consider restricting Hooke to only have 'velocity' as an option for damping? + if (damping[itype][jtype] == VELOCITY){ + damp_normal = normal_coeffs[itype][jtype][1]; + } + else if (damping[itype][jtype] == VISCOELASTIC){ + if (normal[itype][jtype] == HOOKE) a = sqrt(dR); + damp_normal = normal_coeffs[itype][jtype][1]*a*meff; + } + else if (damping[itype][jtype] == TSUJI){ + damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); + } + + Fdamp = -damp_normal*vnnr; + + Fntot = Fne + Fdamp; + + //**************************************** + //Tangential force, including history effects + //**************************************** + + // tangential component + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + wr1 = (radi*omega[i][0] + radj*omega[j][0]); + wr2 = (radi*omega[i][1] + radj*omega[j][1]); + wr3 = (radi*omega[i][2] + radj*omega[j][2]); + + // relative tangential velocities + vtr1 = vt1 - (nz*wr2-ny*wr3); + vtr2 = vt2 - (nx*wr3-nz*wr1); + vtr3 = vt3 - (ny*wr1-nx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + // If any history is needed: + if (use_history){ + touch[jj] = 1; + history = &allhistory[size_history*jj]; + } + + if (normal[itype][jtype] == JKR){ + F_pulloff = 3*M_PI*coh*Reff; + Fcrit = fabs(Fne + 2*F_pulloff); + } + else{ + Fcrit = fabs(Fne); + } + + //------------------------------ + //Tangential forces + //------------------------------ + k_tangential = tangential_coeffs[itype][jtype][0]; + damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; + + if (tangential_history){ + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + + // Rotate and update displacements. + // See e.g. eq. 17 of Luding, Gran. Matter 2008, v10,p235 + if (historyupdate) { + rsht = history[0]*nx + history[1]*ny + history[2]*nz; + if (fabs(rsht) < EPSILON) rsht = 0; + if (rsht > 0){ + scalefac = shrmag/(shrmag - rsht); //if rhst == shrmag, contacting pair has rotated 90 deg. in one step, in which case you deserve a crash! + history[0] -= rsht*nx; + history[1] -= rsht*ny; + history[2] -= rsht*nz; + //Also rescale to preserve magnitude + history[0] *= scalefac; + history[1] *= scalefac; + history[2] *= scalefac; + } + //Update history + history[0] += vtr1*dt; + history[1] += vtr2*dt; + history[2] += vtr3*dt; + } + + // tangential forces = history + tangential velocity damping + fs1 = -k_tangential*history[0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + + // rescale frictional displacements and forces if needed + Fscrit = tangential_coeffs[itype][jtype][2] * Fcrit; + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + if (fs > Fscrit) { + if (shrmag != 0.0) { + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); + history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); + history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); + fs1 *= Fscrit/fs; + fs2 *= Fscrit/fs; + fs3 *= Fscrit/fs; + } else fs1 = fs2 = fs3 = 0.0; + } + } + else{ //Classic pair gran/hooke (no history) + fs = meff*damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; + else Ft = 0.0; + fs1 = -Ft*vtr1; + fs2 = -Ft*vtr2; + fs3 = -Ft*vtr3; + } + + //**************************************** + // Rolling resistance + //**************************************** + + if (roll[itype][jtype] != ROLL_NONE){ + relrot1 = omega[i][0] - omega[j][0]; + relrot2 = omega[i][1] - omega[j][1]; + relrot3 = omega[i][2] - omega[j][2]; + + // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) + // This is different from the Marshall papers, which use the Bagi/Kuhn formulation + // for rolling velocity (see Wang et al for why the latter is wrong) + vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; + vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; + vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; + vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); + if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; + else vrlmaginv = 0.0; + + if (roll_history){ + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; + + // Rolling displacement + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + + if (historyupdate){ + if (fabs(rolldotn) < EPSILON) rolldotn = 0; + if (rolldotn > 0){ //Rotate into tangential plane + scalefac = rollmag/(rollmag - rolldotn); + history[rhist0] -= rolldotn*nx; + history[rhist1] -= rolldotn*ny; + history[rhist2] -= rolldotn*nz; + //Also rescale to preserve magnitude + history[rhist0] *= scalefac; + history[rhist1] *= scalefac; + history[rhist2] *= scalefac; + } + history[rhist0] += vrl1*dt; + history[rhist1] += vrl2*dt; + history[rhist2] += vrl3*dt; + } + + + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][1]; + fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; + fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; + fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; + + // rescale frictional displacements and forces if needed + Frcrit = roll_coeffs[itype][jtype][2] * Fcrit; + + fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); + if (fr > Frcrit) { + if (rollmag != 0.0) { + history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); + history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); + history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); + fr1 *= Frcrit/fr; + fr2 *= Frcrit/fr; + fr3 *= Frcrit/fr; + } else fr1 = fr2 = fr3 = 0.0; + } + } + else{ // + fr = meff*roll_coeffs[itype][jtype][1]*vrlmag; + if (vrlmag != 0.0) fr = MIN(Fne, fr) / vrlmag; + else fr = 0.0; + fr1 = -fr*vrl1; + fr2 = -fr*vrl2; + fr3 = -fr*vrl3; + } + } + + //**************************************** + // Twisting torque, including history effects + //**************************************** + if (twist[itype][jtype] != TWIST_NONE){ + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) + if (twist[itype][jtype] == TWIST_MARSHALL){ + k_twist = 0.5*k_tangential*a*a;; //eq 32 + damp_twist = 0.5*damp_tangential*a*a; + mu_twist = TWOTHIRDS*a; + } + else{ + k_twist = twist_coeffs[itype][jtype][0]; + damp_twist = twist_coeffs[itype][jtype][1]; + mu_twist = twist_coeffs[itype][jtype][2]; + } + if (twist[itype][jtype] > 1){ + if (historyupdate){ + history[twist_history_index] += magtwist*dt; + } + magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) + signtwist = (magtwist > 0) - (magtwist < 0); + Mtcrit = TWOTHIRDS*a*Fscrit;//critical torque (eq 44) + if (fabs(magtortwist) > Mtcrit) { + history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); + magtortwist = -Mtcrit * signtwist; //eq 34 + } + } + else{ + if (magtwist > 0) magtortwist = -damp_twist*magtwist; + else magtortwist = 0; + } + } + // Apply forces & torques + + fx = nx*Fntot + fs1; + fy = ny*Fntot + fs2; + fz = nz*Fntot + fs3; + + f[i][0] += fx; + f[i][1] += fy; + f[i][2] += fz; + + tor1 = ny*fs3 - nz*fs2; + tor2 = nz*fs1 - nx*fs3; + tor3 = nx*fs2 - ny*fs1; + + torque[i][0] -= radi*tor1; + torque[i][1] -= radi*tor2; + torque[i][2] -= radi*tor3; + + if (twist[itype][jtype] != TWIST_NONE){ + tortwist1 = magtortwist * nx; + tortwist2 = magtortwist * ny; + tortwist3 = magtortwist * nz; + + torque[i][0] += tortwist1; + torque[i][1] += tortwist2; + torque[i][2] += tortwist3; + } + + if (roll[itype][jtype] != ROLL_NONE){ + torroll1 = Reff*(ny*fr3 - nz*fr2); //n cross fr + torroll2 = Reff*(nz*fr1 - nx*fr3); + torroll3 = Reff*(nx*fr2 - ny*fr1); + + torque[i][0] += torroll1; + torque[i][1] += torroll2; + torque[i][2] += torroll3; + } + + if (force->newton_pair || j < nlocal) { + f[j][0] -= fx; + f[j][1] -= fy; + f[j][2] -= fz; + + torque[j][0] -= radj*tor1; + torque[j][1] -= radj*tor2; + torque[j][2] -= radj*tor3; + + if (twist[itype][jtype] != TWIST_NONE){ + torque[j][0] -= tortwist1; + torque[j][1] -= tortwist2; + torque[j][2] -= tortwist3; + } + if (roll[itype][jtype] != ROLL_NONE){ + torque[j][0] -= torroll1; + torque[j][1] -= torroll2; + torque[j][2] -= torroll3; + } + } + if (evflag) ev_tally_xyz(i,j,nlocal,0, + 0.0,0.0,fx,fy,fz,delx,dely,delz); + } + } + } +} + + +/* ---------------------------------------------------------------------- +allocate all arrays +------------------------------------------------------------------------- */ + +void PairGranularMulti::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(cut,n+1,n+1,"pair:cut"); + memory->create(normal_coeffs,n+1,n+1,4,"pair:normal_coeffs"); + memory->create(tangential_coeffs,n+1,n+1,3,"pair:tangential_coeffs"); + memory->create(roll_coeffs,n+1,n+1,3,"pair:roll_coeffs"); + memory->create(twist_coeffs,n+1,n+1,3,"pair:twist_coeffs"); + + memory->create(normal,n+1,n+1,"pair:normal"); + memory->create(damping,n+1,n+1,"pair:damping"); + memory->create(tangential,n+1,n+1,"pair:tangential"); + memory->create(roll,n+1,n+1,"pair:roll"); + memory->create(twist,n+1,n+1,"pair:twist"); + + onerad_dynamic = new double[n+1]; + onerad_frozen = new double[n+1]; + maxrad_dynamic = new double[n+1]; + maxrad_frozen = new double[n+1]; +} + +/* ---------------------------------------------------------------------- + global settings +------------------------------------------------------------------------- */ + +void PairGranularMulti::settings(int narg, char **arg) +{ + if (narg < 5) error->all(FLERR,"Illegal pair_style command"); + + int iarg = 0; + + //Some defaults + normal_global = HERTZ; + damping_global = VISCOELASTIC; + tangential_global = TANGENTIAL_MINDLIN; + roll_global = ROLL_NONE; + twist_global = TWIST_NONE; + + tangential_history = 1; + roll_history = twist_history = 0; + + int normal_set, tangential_set; + normal_set = tangential_set = 0; + + while (iarg < narg){ + if (strcmp(arg[iarg], "hooke") == 0){ + if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hooke option"); + normal_global = HOOKE; + normal_set = 1; + memory->create(normal_coeffs_global, 2, "pair:normal_coeffs_global"); + normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping + iarg += 3; + } + else if (strcmp(arg[iarg], "hertz") == 0){ + int num_coeffs = 2; + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); + normal_global = HERTZ; + normal_set = 1; + memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); + normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping + iarg += num_coeffs+1; + } + else if (strcmp(arg[iarg], "hertz/material") == 0){ + int num_coeffs = 3; + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); + normal_global = HERTZ_MATERIAL; + normal_set = 1; + memory->create(normal_coeffs_global, num_coeffs, "pair:normal_coeffs_global"); + normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E (Young's modulus) + normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G (shear modulus) + iarg += num_coeffs+1; + } + else if (strcmp(arg[iarg], "dmt") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for Hertz option"); + normal_global = DMT; + normal_set = 1; + memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); + normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //4/3 E + normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G + normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion + iarg += 5; + } + else if (strcmp(arg[iarg], "jkr") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for JKR option"); + beyond_contact = 1; + normal_global = JKR; + normal_set = 1; + memory->create(normal_coeffs_global, 4, "pair:normal_coeffs_global"); + normal_coeffs_global[0] = force->numeric(FLERR,arg[iarg+1]); //E + normal_coeffs_global[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_global[2] = force->numeric(FLERR,arg[iarg+3]); //G + normal_coeffs_global[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + iarg += 5; + } + else if (strcmp(arg[iarg], "damp_velocity") == 0){ + damping_global = VELOCITY; + iarg += 1; + } + else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ + damping_global = VISCOELASTIC; + iarg += 1; + } + else if (strcmp(arg[iarg], "damp_tsuji") == 0){ + damping_global = TSUJI; + iarg += 1; + } + else if (strstr(arg[iarg], "tangential") != NULL){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for tangential model"); + if (strstr(arg[iarg+1], "nohistory") != NULL){ + tangential_global = TANGENTIAL_NOHISTORY; + tangential_set = 1; + } + else if (strstr(arg[iarg+1], "mindlin") != NULL){ + tangential_global = TANGENTIAL_MINDLIN; + tangential_history = 1; + tangential_set = 1; + } + else{ + error->all(FLERR, "Illegal pair_style command, unrecognized sliding friction model"); + } + memory->create(tangential_coeffs_global, 3, "pair:tangential_coeffs_global"); + tangential_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kt + tangential_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + tangential_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + else if (strstr(arg[iarg], "roll") != NULL){ + if (strstr(arg[iarg+1], "none") != NULL){ + roll_global = ROLL_NONE; + iarg += 2; + } + else{ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for rolling model"); + if (strstr(arg[iarg+1], "nohistory") != NULL){ + roll_global = ROLL_NOHISTORY; + } + else if (strstr(arg[iarg+1], "sds") != NULL){ + roll_global = ROLL_SDS; + roll_history = 1; + } + else{ + error->all(FLERR, "Illegal pair_style command, unrecognized rolling friction model"); + } + memory->create(roll_coeffs_global, 3, "pair:roll_coeffs_global"); + roll_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //kR + roll_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR + roll_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + } + else if (strstr(arg[iarg], "twist") != NULL){ + if (strstr(arg[iarg+1], "none") != NULL){ + twist_global = TWIST_NONE; + iarg += 2; + } + else if (strstr(arg[iarg+1], "marshall") != NULL){ + twist_global = TWIST_MARSHALL; + twist_history = 1; + iarg += 2; + } + else{ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_style command, not enough parameters provided for twist model"); + memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); //To be filled later + if (strstr(arg[iarg+1], "nohistory") != NULL){ + twist_global = TWIST_NOHISTORY; + } + else if (strstr(arg[iarg+1], "sds") != NULL){ + twist_global = TWIST_SDS; + twist_history = 1; + } + else{ + error->all(FLERR, "Illegal pair_style command, unrecognized twisting friction model"); + } + memory->create(twist_coeffs_global, 3, "pair:twist_coeffs_global"); + twist_coeffs_global[0] = force->numeric(FLERR,arg[iarg+2]); //ktwist + twist_coeffs_global[1] = force->numeric(FLERR,arg[iarg+3]); //gammatwist + twist_coeffs_global[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + } + else error->all(FLERR, "Illegal pair_style granular command"); + } + + //Set all i-i entries, which may be replaced by pair coeff commands + //It may also make sense to consider removing all of the above, and only + // having the option for pair_coeff to set the parameters, similar to most LAMMPS pair styles + // The reason for the current setup is to remain true to existing pair gran/hooke etc. syntax, + // where coeffs are set in the pair_style command, and a pair_coeff * * command is issued. + + //Other option is to have two pair styles, e.g. pair granular and pair granular/multi, + // where granular/multi allows per-type coefficients, pair granular does not (this would also + // allow minor speed-up by templating pair granular) + allocate(); + double damp; + for (int i = 1; i <= atom->ntypes; i++){ + normal[i][i] = normal_global; + damping[i][i] = damping_global; + tangential[i][i] = tangential_global; + roll[i][i] = roll_global; + twist[i][i] = twist_global; + + if (normal_set){ + if (damping_global == TSUJI){ + double cor = normal_coeffs_global[1]; + damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ + 27.467*pow(cor,4)-18.022*pow(cor,5)+ + 4.8218*pow(cor,6); + } + else damp = normal_coeffs_global[1]; + normal_coeffs[i][i][0] = normal_coeffs_global[0]; + normal_coeffs[i][i][1] = damp; + if (normal[i][i] != HOOKE && normal[i][i] != HERTZ){ + normal_coeffs[i][i][2] = normal_coeffs_global[2]; + } + if ((normal_global == JKR) || (normal_global == DMT)) + normal_coeffs[i][i][3] = normal_coeffs_global[3]; + } + if(tangential_set){ + tangential[i][i] = tangential_global; + for (int k = 0; k < 3; k++) + tangential_coeffs[i][i][k] = tangential_coeffs_global[k]; + } + roll[i][i] = roll_global; + if (roll_global != ROLL_NONE) + for (int k = 0; k < 3; k++) + roll_coeffs[i][i][k] = roll_coeffs_global[k]; + + twist[i][i] = twist_global; + if (twist_global != TWIST_NONE && twist_global != TWIST_MARSHALL) + for (int k = 0; k < 3; k++) + twist_coeffs[i][i][k] = twist_coeffs_global[k]; + + if (normal_set && tangential_set) setflag[i][i] = 1; + } +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairGranularMulti::coeff(int narg, char **arg) +{ + int normal_local, damping_local, tangential_local, roll_local, twist_local; + double *normal_coeffs_local; + double *tangential_coeffs_local; + double *roll_coeffs_local; + double *twist_coeffs_local; + + normal_coeffs_local = new double[4]; + tangential_coeffs_local = new double[4]; + roll_coeffs_local = new double[4]; + twist_coeffs_local = new double[4]; + + if (narg < 2) + error->all(FLERR,"Incorrect args for pair coefficients"); + + if (!allocated) allocate(); + + int ilo,ihi,jlo,jhi; + force->bounds(FLERR,arg[0],atom->ntypes,ilo,ihi); + force->bounds(FLERR,arg[1],atom->ntypes,jlo,jhi); + + normal_local = tangential_local = roll_local = twist_local = -1; + damping_local = -1; + + int iarg = 2; + while (iarg < narg){ + if (strcmp(arg[iarg], "hooke") == 0){ + if (iarg + 2 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hooke option"); + normal_local = HOOKE; + normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + iarg += 3; + } + else if (strcmp(arg[iarg], "hertz") == 0){ + int num_coeffs = 2; + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); + normal_local = HERTZ; + normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + iarg += num_coeffs+1; + } + else if (strcmp(arg[iarg], "hertz/material") == 0){ + int num_coeffs = 3; + if (iarg + num_coeffs >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); + normal_local = HERTZ; + normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E + normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G + iarg += num_coeffs+1; + } + else if (strcmp(arg[iarg], "dmt") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for Hertz option"); + normal_local = DMT; + normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E + normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G + normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+3]); //cohesion + iarg += 5; + } + else if (strcmp(arg[iarg], "jkr") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for JKR option"); + beyond_contact = 1; + normal_local = JKR; + normal_coeffs_local[0] = force->numeric(FLERR,arg[iarg+1]); //E + normal_coeffs_local[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_local[2] = force->numeric(FLERR,arg[iarg+3]); //G + normal_coeffs_local[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion + iarg += 5; + } + else if (strcmp(arg[iarg], "damp_velocity") == 0){ + damping_local = VELOCITY; + iarg += 1; + } + else if (strcmp(arg[iarg], "damp_viscoelastic") == 0){ + damping_local = VISCOELASTIC; + iarg += 1; + } + else if (strcmp(arg[iarg], "damp_tsuji") == 0){ + damping_local = TSUJI; + iarg += 1; + } + else if (strstr(arg[iarg], "tangential") != NULL){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + if (strstr(arg[iarg+1], "nohistory") != NULL){ + tangential_local = TANGENTIAL_NOHISTORY; + } + else if (strstr(arg[iarg+1], "mindlin") != NULL){ + tangential_local = TANGENTIAL_MINDLIN; + tangential_history = 1; + } + else{ + error->all(FLERR, "Illegal pair_coeff command, tangential model not recognized"); + } + tangential_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt + tangential_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + tangential_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + else if (strstr(arg[iarg], "rolling") != NULL){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (strstr(arg[iarg+1], "none") != NULL){ + roll_local = ROLL_NONE; + iarg += 2; + } + else{ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); + if (strstr(arg[iarg+1], "nohistory") != NULL){ + roll_local = ROLL_NOHISTORY; + } + else if (strstr(arg[iarg+1], "sds") != NULL){ + roll_local = ROLL_SDS; + roll_history = 1; + } + else{ + error->all(FLERR, "Illegal pair_coeff command, rolling friction model not recognized"); + } + roll_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt + roll_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + roll_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + } + else if (strstr(arg[iarg], "twist") != NULL){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (strstr(arg[iarg+1], "none") != NULL){ + twist_local = TWIST_NONE; + iarg += 2; + } + else if (strstr(arg[iarg+1], "marshall") != NULL){ + twist_local = TWIST_MARSHALL; + twist_history = 1; + iarg += 2; + } + else{ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); + if (strstr(arg[iarg+1], "nohistory") != NULL){ + twist_local = TWIST_NOHISTORY; + } + else if (strstr(arg[iarg+1], "sds") != NULL){ + twist_local = TWIST_SDS; + twist_history = 1; + } + else{ + error->all(FLERR, "Illegal pair_coeff command, twisting friction model not recognized"); + } + twist_coeffs_local[0] = force->numeric(FLERR,arg[iarg+2]); //kt + twist_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + twist_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. + iarg += 5; + } + } + else error->all(FLERR, "Illegal pair coeff command"); + } + + int count = 0; + double damp; + if (damping_local >= 0){ + if (normal_local == -1) + error->all(FLERR, "Illegal pair_coeff command, must specify normal model when setting damping model"); + if (damping_local == TSUJI){ + double cor; + cor = normal_coeffs_local[1]; + damp = 1.2728-4.2783*cor+11.087*pow(cor,2)-22.348*pow(cor,3)+ + 27.467*pow(cor,4)-18.022*pow(cor,5)+ + 4.8218*pow(cor,6); + } + else damp = normal_coeffs_local[1]; + } + + for (int i = ilo; i <= ihi; i++) { + for (int j = MAX(jlo,i); j <= jhi; j++) { + if (normal_local >= 0){ + normal[i][j] = normal_local; + normal_coeffs[i][j][0] = normal_coeffs_local[0]; + if (damping_local == -1){ + damp = normal_coeffs_global[1]; + } + normal_coeffs[i][j][1] = damp; + if (normal_local != HERTZ && normal_local != HOOKE) normal_coeffs[i][j][2] = normal_coeffs_local[2]; + if ((normal_local == JKR) || (normal_local == DMT)) + normal_coeffs[i][j][3] = normal_coeffs_local[3]; + } + if (damping_local >= 0){ + damping[i][j] = damping_local; + } + if (tangential_local >= 0){ + tangential[i][j] = tangential_local; + for (int k = 0; k < 3; k++) + tangential_coeffs[i][j][k] = tangential_coeffs_local[k]; + } + if (roll_local >= 0){ + roll[i][j] = roll_local; + if (roll_local != ROLL_NONE) + for (int k = 0; k < 3; k++) + roll_coeffs[i][j][k] = roll_coeffs_local[k]; + } + if (twist_local >= 0){ + twist[i][j] = twist_local; + if (twist_local != TWIST_NONE && twist_local != TWIST_MARSHALL) + for (int k = 0; k < 3; k++) + twist_coeffs[i][j][k] = twist_coeffs_local[k]; + } + + if (normal_local >= 0 && tangential_local >= 0) setflag[i][j] = 1; + + count++; + } + } + + delete[] normal_coeffs_local; + delete[] tangential_coeffs_local; + delete[] roll_coeffs_local; + delete[] twist_coeffs_local; + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairGranularMulti::init_style() +{ + int i; + + // error and warning checks + + if (!atom->radius_flag || !atom->rmass_flag) + error->all(FLERR,"Pair granular requires atom attributes radius, rmass"); + if (comm->ghost_velocity == 0) + error->all(FLERR,"Pair granular requires ghost atoms store velocity"); + + // Determine whether we need a granular neigh list, how large it needs to be + use_history = tangential_history || roll_history || twist_history; + + //For JKR, will need fix/neigh/history to keep track of touch arrays + for (int i = 1; i <= atom->ntypes; i++) + for (int j = 1; j <= atom->ntypes; j++) + if (normal[i][j] == JKR) use_history = 1; + + size_history = 3*tangential_history + 3*roll_history + twist_history; + + //Determine location of tangential/roll/twist histories in array + if (roll_history){ + if (tangential_history) roll_history_index = 3; + else roll_history_index = 0; + } + if (twist_history){ + if (tangential_history){ + if (roll_history) twist_history_index = 6; + else twist_history_index = 3; + } + else{ + if (roll_history) twist_history_index = 3; + else twist_history_index = 0; + } + } + + int irequest = neighbor->request(this,instance_me); + neighbor->requests[irequest]->size = 1; + if (use_history) neighbor->requests[irequest]->history = 1; + + dt = update->dt; + + // if history is stored: + // if first init, create Fix needed for storing history + + if (use_history && fix_history == NULL) { + char dnumstr[16]; + sprintf(dnumstr,"%d",size_history); + char **fixarg = new char*[4]; + fixarg[0] = (char *) "NEIGH_HISTORY"; + fixarg[1] = (char *) "all"; + fixarg[2] = (char *) "NEIGH_HISTORY"; + fixarg[3] = dnumstr; + modify->add_fix(4,fixarg,1); + delete [] fixarg; + fix_history = (FixNeighHistory *) modify->fix[modify->nfix-1]; + fix_history->pair = this; + } + + // check for FixFreeze and set freeze_group_bit + + for (i = 0; i < modify->nfix; i++) + if (strcmp(modify->fix[i]->style,"freeze") == 0) break; + if (i < modify->nfix) freeze_group_bit = modify->fix[i]->groupbit; + else freeze_group_bit = 0; + + // check for FixRigid so can extract rigid body masses + + fix_rigid = NULL; + for (i = 0; i < modify->nfix; i++) + if (modify->fix[i]->rigid_flag) break; + if (i < modify->nfix) fix_rigid = modify->fix[i]; + + // check for FixPour and FixDeposit so can extract particle radii + + int ipour; + for (ipour = 0; ipour < modify->nfix; ipour++) + if (strcmp(modify->fix[ipour]->style,"pour") == 0) break; + if (ipour == modify->nfix) ipour = -1; + + int idep; + for (idep = 0; idep < modify->nfix; idep++) + if (strcmp(modify->fix[idep]->style,"deposit") == 0) break; + if (idep == modify->nfix) idep = -1; + + // set maxrad_dynamic and maxrad_frozen for each type + // include future FixPour and FixDeposit particles as dynamic + + int itype; + for (i = 1; i <= atom->ntypes; i++) { + onerad_dynamic[i] = onerad_frozen[i] = 0.0; + if (ipour >= 0) { + itype = i; + double radmax = *((double *) modify->fix[ipour]->extract("radius",itype)); + if (normal[itype][itype] == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); + onerad_dynamic[i] = radmax; + } + if (idep >= 0) { + itype = i; + double radmax = *((double *) modify->fix[idep]->extract("radius",itype)); + if (normal[itype][itype] == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); + onerad_dynamic[i] = radmax; + } + } + + double *radius = atom->radius; + int *mask = atom->mask; + int *type = atom->type; + int nlocal = atom->nlocal; + + for (i = 0; i < nlocal; i++){ + double radius_cut = radius[i]; + if (normal[type[i]][type[i]] == JKR){ + radius_cut = radius[i] - 0.5*pulloff_distance(radius[i], type[i]); + } + if (mask[i] & freeze_group_bit){ + onerad_frozen[type[i]] = MAX(onerad_frozen[type[i]],radius_cut); + } + else{ + onerad_dynamic[type[i]] = MAX(onerad_dynamic[type[i]],radius_cut); + } + } + + MPI_Allreduce(&onerad_dynamic[1],&maxrad_dynamic[1],atom->ntypes, + MPI_DOUBLE,MPI_MAX,world); + MPI_Allreduce(&onerad_frozen[1],&maxrad_frozen[1],atom->ntypes, + MPI_DOUBLE,MPI_MAX,world); + + // set fix which stores history info + + if (size_history > 0){ + int ifix = modify->find_fix("NEIGH_HISTORY"); + if (ifix < 0) error->all(FLERR,"Could not find pair fix neigh history ID"); + fix_history = (FixNeighHistory *) modify->fix[ifix]; + } +} + +/* ---------------------------------------------------------------------- + init for one type pair i,j and corresponding j,i +------------------------------------------------------------------------- */ + +double PairGranularMulti::init_one(int i, int j) +{ + double cutoff; + if (setflag[i][j] == 0) { + if ((normal[i][i] != normal[j][j]) || + (damping[i][i] != damping[j][j]) || + (tangential[i][i] != tangential[j][j]) || + (roll[i][i] != roll[j][j]) || + (twist[i][i] != twist[j][j])){ + + char str[512]; + sprintf(str,"Granular pair style functional forms are different, cannot mix coefficients for types %d and %d. \nThis combination must be set explicitly via pair_coeff command.",i,j); + error->one(FLERR,str); + } + + if (normal[i][j] != HOOKE && normal[i][j] != HERTZ){ + normal_coeffs[i][j][0] = mix_stiffnessE(normal_coeffs[i][i][0], normal_coeffs[j][j][0], + normal_coeffs[i][i][2], normal_coeffs[j][j][2]); + normal_coeffs[i][j][2] = mix_stiffnessG(normal_coeffs[i][i][0], normal_coeffs[j][j][0], + normal_coeffs[i][i][2], normal_coeffs[j][j][2]); + } + else{ + normal_coeffs[i][j][0] = mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); + } + + normal_coeffs[i][j][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); + if ((normal[i][i] == JKR) || (normal[i][i] == DMT)) + normal_coeffs[i][j][3] = mix_geom(normal_coeffs[i][i][3], normal_coeffs[j][j][3]); + + for (int k = 0; k < 3; k++) + tangential_coeffs[i][j][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); + + + if (roll[i][i] != ROLL_NONE){ + for (int k = 0; k < 3; k++) + roll_coeffs[i][j][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); + } + + if (twist[i][i] != TWIST_NONE && twist[i][i] != TWIST_MARSHALL){ + for (int k = 0; k < 3; k++) + twist_coeffs[i][j][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]); + } + } + + // It is possible that cut[i][j] at this point is still 0.0. This can happen when + // there is a future fix_pour after the current run. A cut[i][j] = 0.0 creates + // problems because neighbor.cpp uses min(cut[i][j]) to decide on the bin size + // To avoid this issue, for cases involving cut[i][j] = 0.0 (possible only + // if there is no current information about radius/cutoff of type i and j). + // we assign cutoff = max(cut[i][j]) for i,j such that cut[i][j] > 0.0. + + if (((maxrad_dynamic[i] > 0.0) && (maxrad_dynamic[j] > 0.0)) || + ((maxrad_dynamic[i] > 0.0) && (maxrad_frozen[j] > 0.0)) || + ((maxrad_frozen[i] > 0.0) && (maxrad_dynamic[j] > 0.0))) { // radius info about both i and j exist + cutoff = maxrad_dynamic[i]+maxrad_dynamic[j]; + cutoff = MAX(cutoff,maxrad_frozen[i]+maxrad_dynamic[j]); + cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]); + } + else { // radius info about either i or j does not exist (i.e. not present and not about to get poured; set to largest value to not interfere with neighbor list) + double cutmax = 0.0; + for (int k = 1; k <= atom->ntypes; k++) { + cutmax = MAX(cutmax,2.0*maxrad_dynamic[k]); + cutmax = MAX(cutmax,2.0*maxrad_frozen[k]); + } + cutoff = cutmax; + } + return cutoff; +} + + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file + ------------------------------------------------------------------------- */ + +void PairGranularMulti::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(&normal[i][j],sizeof(int),1,fp); + fwrite(&damping[i][j],sizeof(int),1,fp); + fwrite(&tangential[i][j],sizeof(int),1,fp); + fwrite(&roll[i][j],sizeof(int),1,fp); + fwrite(&twist[i][j],sizeof(int),1,fp); + fwrite(&normal_coeffs[i][j],sizeof(double),4,fp); + fwrite(&tangential_coeffs[i][j],sizeof(double),3,fp); + fwrite(&roll_coeffs[i][j],sizeof(double),3,fp); + fwrite(&twist_coeffs[i][j],sizeof(double),3,fp); + fwrite(&cut[i][j],sizeof(double),1,fp); + } + } + } +} + +/* ---------------------------------------------------------------------- + proc 0 reads from restart file, bcasts + ------------------------------------------------------------------------- */ + +void PairGranularMulti::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(&normal[i][j],sizeof(int),1,fp); + fread(&damping[i][j],sizeof(int),1,fp); + fread(&tangential[i][j],sizeof(int),1,fp); + fread(&roll[i][j],sizeof(int),1,fp); + fread(&twist[i][j],sizeof(int),1,fp); + fread(&normal_coeffs[i][j],sizeof(double),4,fp); + fread(&tangential_coeffs[i][j],sizeof(double),3,fp); + fread(&roll_coeffs[i][j],sizeof(double),3,fp); + fread(&twist_coeffs[i][j],sizeof(double),3,fp); + fread(&cut[i][j],sizeof(double),1,fp); + } + MPI_Bcast(&normal[i][j],1,MPI_INT,0,world); + MPI_Bcast(&damping[i][j],1,MPI_INT,0,world); + MPI_Bcast(&tangential[i][j],1,MPI_INT,0,world); + MPI_Bcast(&roll[i][j],1,MPI_INT,0,world); + MPI_Bcast(&twist[i][j],1,MPI_INT,0,world); + MPI_Bcast(&normal_coeffs[i][j],4,MPI_DOUBLE,0,world); + MPI_Bcast(&tangential_coeffs[i][j],3,MPI_DOUBLE,0,world); + MPI_Bcast(&roll_coeffs[i][j],3,MPI_DOUBLE,0,world); + MPI_Bcast(&twist_coeffs[i][j],3,MPI_DOUBLE,0,world); + MPI_Bcast(&cut[i][j],1,MPI_DOUBLE,0,world); + } + } + } +} + + +/* ---------------------------------------------------------------------- */ + +void PairGranularMulti::reset_dt() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +double PairGranularMulti::single(int i, int j, int itype, int jtype, + double rsq, double factor_coul, double factor_lj, double &fforce) +{ + double radi,radj,radsum; + double r,rinv,rsqinv,delx,dely,delz, nx, ny, nz, Reff; + double dR, dR2; + double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3,wr1,wr2,wr3; + double vtr1,vtr2,vtr3,vrel; + double mi,mj,meff,damp,ccel,tor1,tor2,tor3; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; + + double knfac, damp_normal; + double k_tangential, damp_tangential; + double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; + double fs, fs1, fs2, fs3; + + //For JKR + double R2, coh, F_pulloff, delta_pulloff, dist_pulloff, a, a2, E; + double delta, t0, t1, t2, t3, t4, t5, t6; + double sqrt1, sqrt2, sqrt3, sqrt4; + + + //Rolling + double k_roll, damp_roll; + double roll1, roll2, roll3, torroll1, torroll2, torroll3; + double rollmag, rolldotn, scalefac; + double fr, fr1, fr2, fr3; + + //Twisting + double k_twist, damp_twist, mu_twist; + double signtwist, magtwist, magtortwist, Mtcrit; + double tortwist1, tortwist2, tortwist3; + + double shrmag,rsht; + int jnum; + int *ilist,*jlist,*numneigh,**firstneigh; + int *touch,**firsttouch; + double *history,*allhistory,**firsthistory; + + double *radius = atom->radius; + radi = radius[i]; + radj = radius[j]; + radsum = radi + radj; + Reff = radi*radj/(radi+radj); + + bool touchflag; + if (normal[itype][jtype] == JKR){ + R2 = Reff*Reff; + coh = normal_coeffs[itype][jtype][3]; + a = cbrt(9.0*M_PI*coh*R2/(4*E)); + delta_pulloff = a*a/Reff - 2*sqrt(M_PI*coh*a/E); + dist_pulloff = radsum+delta_pulloff; + touchflag = (rsq <= dist_pulloff*dist_pulloff); + } + else{ + touchflag = (rsq <= radsum*radsum); + } + + if (touchflag){ + fforce = 0.0; + for (int m = 0; m < single_extra; m++) svector[m] = 0.0; + return 0.0; + } + + double **x = atom->x; + delx = x[i][0] - x[j][0]; + dely = x[i][1] - x[j][1]; + delz = x[i][2] - x[j][2]; + r = sqrt(rsq); + rinv = 1.0/r; + + nx = delx*rinv; + ny = dely*rinv; + nz = delz*rinv; + + // relative translational velocity + + double **v = atom->v; + vr1 = v[i][0] - v[j][0]; + vr2 = v[i][1] - v[j][1]; + vr3 = v[i][2] - v[j][2]; + + // normal component + + vnnr = vr1*nx + vr2*ny + vr3*nz; + vn1 = nx*vnnr; + vn2 = ny*vnnr; + vn3 = nz*vnnr; + + double *rmass = atom->rmass; + int *mask = atom->mask; + mi = rmass[i]; + mj = rmass[j]; + if (fix_rigid) { + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; + } + + meff = mi*mj / (mi+mj); + if (mask[i] & freeze_group_bit) meff = mj; + if (mask[j] & freeze_group_bit) meff = mi; + + delta = radsum - r; + dR = delta*Reff; + + // tangential component + + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + + double **omega = atom->omega; + wr1 = (radi*omega[i][0] + radj*omega[j][0]); + wr2 = (radi*omega[i][1] + radj*omega[j][1]); + wr3 = (radi*omega[i][2] + radj*omega[j][2]); + + // meff = effective mass of pair of particles + // if I or J part of rigid body, use body mass + // if I or J is frozen, meff is other particle + + int *type = atom->type; + + mi = rmass[i]; + mj = rmass[j]; + if (fix_rigid) { + // NOTE: ensure mass_rigid is current for owned+ghost atoms? + if (mass_rigid[i] > 0.0) mi = mass_rigid[i]; + if (mass_rigid[j] > 0.0) mj = mass_rigid[j]; + } + + meff = mi*mj / (mi+mj); + if (mask[i] & freeze_group_bit) meff = mj; + if (mask[j] & freeze_group_bit) meff = mi; + + delta = radsum - r; + dR = delta*Reff; + if (normal[itype][jtype] == JKR){ + dR2 = dR*dR; + t0 = coh*coh*R2*R2*E; + t1 = PI27SQ*t0; + t2 = 8*dR*dR2*E*E*E; + t3 = 4*dR2*E; + sqrt1 = MAX(0, t0*(t1+2*t2)); //In case of sqrt(0) < 0 due to precision issues + t4 = cbrt(t1+t2+THREEROOT3*M_PI*sqrt(sqrt1)); + t5 = t3/t4 + t4/E; + sqrt2 = MAX(0, 2*dR + t5); + t6 = sqrt(sqrt2); + sqrt3 = MAX(0, 4*dR - t5 + SIXROOT6*coh*M_PI*R2/(E*t6)); + a = INVROOT6*(t6 + sqrt(sqrt3)); + a2 = a*a; + knfac = FOURTHIRDS*E*a; + Fne = knfac*a2/Reff - TWOPI*a2*sqrt(4*coh*E/(M_PI*a)); + } + else{ + knfac = E; + Fne = knfac*delta; + if (normal[itype][jtype] != HOOKE) + a = sqrt(dR); + Fne *= a; + if (normal[itype][jtype] == DMT) + Fne -= 4*MY_PI*normal_coeffs[itype][jtype][3]*Reff; + } + + //Consider restricting Hooke to only have 'velocity' as an option for damping? + if (damping[itype][jtype] == VELOCITY){ + damp_normal = normal_coeffs[itype][jtype][1]; + } + else if (damping[itype][jtype] == VISCOELASTIC){ + if (normal[itype][jtype] == HOOKE) a = sqrt(dR); + damp_normal = normal_coeffs[itype][jtype][1]*a*meff; + } + else if (damping[itype][jtype] == TSUJI){ + damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); + } + + Fdamp = -damp_normal*vnnr; + + Fntot = Fne + Fdamp; + + jnum = list->numneigh[i]; + jlist = list->firstneigh[i]; + + if (use_history){ + allhistory = fix_history->firstvalue[i]; + for (int jj = 0; jj < jnum; jj++) { + neighprev++; + if (neighprev >= jnum) neighprev = 0; + if (jlist[neighprev] == j) break; + } + history = &allhistory[size_history*neighprev]; + } + + //**************************************** + //Tangential force, including history effects + //**************************************** + + // tangential component + vt1 = vr1 - vn1; + vt2 = vr2 - vn2; + vt3 = vr3 - vn3; + + // relative rotational velocity + wr1 = (radi*omega[i][0] + radj*omega[j][0]); + wr2 = (radi*omega[i][1] + radj*omega[j][1]); + wr3 = (radi*omega[i][2] + radj*omega[j][2]); + + // relative tangential velocities + vtr1 = vt1 - (nz*wr2-ny*wr3); + vtr2 = vt2 - (nx*wr3-nz*wr1); + vtr3 = vt3 - (ny*wr1-nx*wr2); + vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; + vrel = sqrt(vrel); + + if (normal[itype][jtype] == JKR){ + F_pulloff = 3*M_PI*coh*Reff; + Fcrit = fabs(Fne + 2*F_pulloff); + } + else{ + Fcrit = fabs(Fne); + } + + //------------------------------ + //Tangential forces + //------------------------------ + k_tangential = tangential_coeffs[itype][jtype][0]; + if (normal[itype][jtype] != HOOKE){ + k_tangential *= a; + } + damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; + + if (tangential_history){ + shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + + history[2]*history[2]); + + // tangential forces = history + tangential velocity damping + fs1 = -k_tangential*history[0] - damp_tangential*vtr1; + fs2 = -k_tangential*history[1] - damp_tangential*vtr2; + fs3 = -k_tangential*history[2] - damp_tangential*vtr3; + + // rescale frictional displacements and forces if needed + Fscrit = tangential_coeffs[itype][jtype][2] * Fcrit; + fs = sqrt(fs1*fs1 + fs2*fs2 + fs3*fs3); + if (fs > Fscrit) { + if (shrmag != 0.0) { + history[0] = -1.0/k_tangential*(Fscrit*fs1/fs + damp_tangential*vtr1); + history[1] = -1.0/k_tangential*(Fscrit*fs2/fs + damp_tangential*vtr2); + history[2] = -1.0/k_tangential*(Fscrit*fs3/fs + damp_tangential*vtr3); + fs1 *= Fscrit/fs; + fs2 *= Fscrit/fs; + fs3 *= Fscrit/fs; + } else fs1 = fs2 = fs3 = 0.0; + } + } + else{ //Classic pair gran/hooke (no history) + fs = meff*damp_tangential*vrel; + if (vrel != 0.0) Ft = MIN(Fne,fs) / vrel; + else Ft = 0.0; + fs1 = -Ft*vtr1; + fs2 = -Ft*vtr2; + fs3 = -Ft*vtr3; + } + + //**************************************** + // Rolling resistance + //**************************************** + + if (roll[itype][jtype] != ROLL_NONE){ + relrot1 = omega[i][0] - omega[j][0]; + relrot2 = omega[i][1] - omega[j][1]; + relrot3 = omega[i][2] - omega[j][2]; + + // rolling velocity, see eq. 31 of Wang et al, Particuology v 23, p 49 (2015) + // This is different from the Marshall papers, which use the Bagi/Kuhn formulation + // for rolling velocity (see Wang et al for why the latter is wrong) + vrl1 = Reff*(relrot2*nz - relrot3*ny); //- 0.5*((radj-radi)/radsum)*vtr1; + vrl2 = Reff*(relrot3*nx - relrot1*nz); //- 0.5*((radj-radi)/radsum)*vtr2; + vrl3 = Reff*(relrot1*ny - relrot2*nx); //- 0.5*((radj-radi)/radsum)*vtr3; + vrlmag = sqrt(vrl1*vrl1+vrl2*vrl2+vrl3*vrl3); + if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; + else vrlmaginv = 0.0; + + if (roll_history){ + int rhist0 = roll_history_index; + int rhist1 = rhist0 + 1; + int rhist2 = rhist1 + 1; + + // Rolling displacement + rollmag = sqrt(history[rhist0]*history[rhist0] + + history[rhist1]*history[rhist1] + + history[rhist2]*history[rhist2]); + + rolldotn = history[rhist0]*nx + history[rhist1]*ny + history[rhist2]*nz; + + k_roll = roll_coeffs[itype][jtype][0]; + damp_roll = roll_coeffs[itype][jtype][1]; + fr1 = -k_roll*history[rhist0] - damp_roll*vrl1; + fr2 = -k_roll*history[rhist1] - damp_roll*vrl2; + fr3 = -k_roll*history[rhist2] - damp_roll*vrl3; + + // rescale frictional displacements and forces if needed + Frcrit = roll_coeffs[itype][jtype][2] * Fcrit; + + fr = sqrt(fr1*fr1 + fr2*fr2 + fr3*fr3); + if (fr > Frcrit) { + if (rollmag != 0.0) { + history[rhist0] = -1.0/k_roll*(Frcrit*fr1/fr + damp_roll*vrl1); + history[rhist1] = -1.0/k_roll*(Frcrit*fr2/fr + damp_roll*vrl2); + history[rhist2] = -1.0/k_roll*(Frcrit*fr3/fr + damp_roll*vrl3); + fr1 *= Frcrit/fr; + fr2 *= Frcrit/fr; + fr3 *= Frcrit/fr; + } else fr1 = fr2 = fr3 = 0.0; + } + } + else{ // + fr = meff*roll_coeffs[itype][jtype][1]*vrlmag; + if (vrlmag != 0.0) fr = MIN(Fne, fr) / vrlmag; + else fr = 0.0; + fr1 = -fr*vrl1; + fr2 = -fr*vrl2; + fr3 = -fr*vrl3; + } + } + + //**************************************** + // Twisting torque, including history effects + //**************************************** + if (twist[itype][jtype] != TWIST_NONE){ + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) + if (twist[itype][jtype] == TWIST_MARSHALL){ + k_twist = 0.5*k_tangential*a*a;; //eq 32 + damp_twist = 0.5*damp_tangential*a*a; + mu_twist = TWOTHIRDS*a; + } + else{ + k_twist = twist_coeffs[itype][jtype][0]; + damp_twist = twist_coeffs[itype][jtype][1]; + mu_twist = twist_coeffs[itype][jtype][2]; + } + if (twist_history){ + magtortwist = -k_twist*history[twist_history_index] - damp_twist*magtwist;//M_t torque (eq 30) + signtwist = (magtwist > 0) - (magtwist < 0); + Mtcrit = TWOTHIRDS*a*Fscrit;//critical torque (eq 44) + if (fabs(magtortwist) > Mtcrit) { + history[twist_history_index] = 1.0/k_twist*(Mtcrit*signtwist - damp_twist*magtwist); + magtortwist = -Mtcrit * signtwist; //eq 34 + } + } + else{ + if (magtwist > 0) magtortwist = -damp_twist*magtwist; + else magtortwist = 0; + } + } + + // set single_extra quantities + + svector[0] = fs1; + svector[1] = fs2; + svector[2] = fs3; + svector[3] = fs; + svector[4] = fr1; + svector[5] = fr2; + svector[6] = fr3; + svector[7] = fr; + svector[8] = magtortwist; + return 0.0; +} + +/* ---------------------------------------------------------------------- */ + +int PairGranularMulti::pack_forward_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++] = mass_rigid[j]; + } + return m; +} + +/* ---------------------------------------------------------------------- */ + +void PairGranularMulti::unpack_forward_comm(int n, int first, double *buf) +{ + int i,m,last; + + m = 0; + last = first + n; + for (i = first; i < last; i++) + mass_rigid[i] = buf[m++]; +} + +/* ---------------------------------------------------------------------- + memory usage of local atom-based arrays + ------------------------------------------------------------------------- */ + +double PairGranularMulti::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + mixing of Young's modulus (E) +------------------------------------------------------------------------- */ + +double PairGranularMulti::mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj) +{ + double poisii = Eii/(2.0*Gii) - 1.0; + double poisjj = Ejj/(2.0*Gjj) - 1.0; + return 1/((1-poisii*poisjj)/Eii+(1-poisjj*poisjj)/Ejj); +} + +/* ---------------------------------------------------------------------- + mixing of shear modulus (G) + ------------------------------------------------------------------------- */ + +double PairGranularMulti::mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj) +{ + double poisii = Eii/(2.0*Gii) - 1.0; + double poisjj = Ejj/(2.0*Gjj) - 1.0; + return 1/((2.0 -poisjj)/Gii+(2.0-poisjj)/Gjj); +} + +/* ---------------------------------------------------------------------- + mixing of everything else +------------------------------------------------------------------------- */ + +double PairGranularMulti::mix_geom(double valii, double valjj) +{ + return sqrt(valii*valjj); +} + + +/* ---------------------------------------------------------------------- + Compute pull-off distance (beyond contact) for a given radius and atom type +------------------------------------------------------------------------- */ + +double PairGranularMulti::pulloff_distance(double radius, int itype) +{ + double E, coh, a, delta_pulloff; + coh = normal_coeffs[itype][itype][3]; + E = mix_stiffnessE(normal_coeffs[itype][itype][0], normal_coeffs[itype][itype][0], + normal_coeffs[itype][itype][2], normal_coeffs[itype][itype][2]); + a = cbrt(9*M_PI*coh*radius*radius/(4*E)); + return a*a/radius - 2*sqrt(M_PI*coh*a/E); +} + diff --git a/src/GRANULAR/pair_granular_multi.h b/src/GRANULAR/pair_granular_multi.h new file mode 100644 index 0000000000..e853e564df --- /dev/null +++ b/src/GRANULAR/pair_granular_multi.h @@ -0,0 +1,108 @@ +/* ---------------------------------------------------------- + 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(granular/multi,PairGranularMulti) + +#else + +#ifndef LMP_PAIR_GRANULAR_MULTI_H +#define LMP_PAIR_GRANULAR_MULTI_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairGranularMulti : public Pair { +public: + PairGranularMulti(class LAMMPS *); + virtual ~PairGranularMulti(); + virtual void compute(int, int); + virtual void settings(int, char **); + virtual void coeff(int, char **); + void init_style(); + double init_one(int, int); + void write_restart(FILE *); + void read_restart(FILE *); + void reset_dt(); + virtual double single(int, int, int, int, double, double, double, double &); + int pack_forward_comm(int, int *, double *, int, int *); + void unpack_forward_comm(int, int, double *); + double memory_usage(); + + protected: + double cut_global; + double dt; + int freeze_group_bit; + int use_history; + + int neighprev; + double *onerad_dynamic,*onerad_frozen; + double *maxrad_dynamic,*maxrad_frozen; + double **cut; + + class FixNeighHistory *fix_history; + + // storage of rigid body masses for use in granular interactions + + class Fix *fix_rigid; // ptr to rigid body fix, NULL if none + double *mass_rigid; // rigid mass for owned+ghost atoms + int nmax; // allocated size of mass_rigid + + virtual void allocate(); + +private: + int size_history; + + //Per-type models + int **normal, **damping, **tangential, **roll, **twist; + + int normal_global, damping_global; + int tangential_global, roll_global, twist_global; + + int tangential_history, roll_history, twist_history; + int tangential_history_index; + int roll_history_index; + int twist_history_index; + + double *normal_coeffs_global; + double *tangential_coeffs_global; + double *roll_coeffs_global; + double *twist_coeffs_global; + + double ***normal_coeffs; + double ***tangential_coeffs; + double ***roll_coeffs; + double ***twist_coeffs; + + double mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj); + double mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj); + double mix_geom(double valii, double valjj); + double pulloff_distance(double radius, int itype); +}; + +} + +#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. + + */ diff --git a/src/fix_neigh_history.cpp b/src/fix_neigh_history.cpp index c21b494aa4..e33ebe57dc 100644 --- a/src/fix_neigh_history.cpp +++ b/src/fix_neigh_history.cpp @@ -408,7 +408,8 @@ void FixNeighHistory::pre_exchange_newton() m = npartner[j]++; partner[j][m] = tag[i]; jvalues = &valuepartner[j][dnum*m]; - for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; + if (pair->nondefault_history_transfer) pair->transfer_history(onevalues, jvalues); + else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } } @@ -520,7 +521,8 @@ void FixNeighHistory::pre_exchange_no_newton() m = npartner[j]++; partner[j][m] = tag[i]; jvalues = &valuepartner[j][dnum*m]; - for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; + if (pair->nondefault_history_transfer) pair->transfer_history(onevalues, jvalues); + else for (n = 0; n < dnum; n++) jvalues[n] = -onevalues[n]; } } } @@ -604,7 +606,7 @@ void FixNeighHistory::post_neighbor() for (jj = 0; jj < jnum; jj++) { j = jlist[jj]; - rflag = sbmask(j); + rflag = sbmask(j) | pair->beyond_contact; j &= NEIGHMASK; jlist[jj] = j; diff --git a/src/pair.h b/src/pair.h index 27b6d41eef..0911bff706 100644 --- a/src/pair.h +++ b/src/pair.h @@ -98,6 +98,8 @@ class Pair : protected Pointers { enum{GEOMETRIC,ARITHMETIC,SIXTHPOWER}; // mixing options + int beyond_contact, nondefault_history_transfer; //for granular styles + // KOKKOS host/device flag and data masks ExecutionSpace execution_space; @@ -180,6 +182,7 @@ class Pair : protected Pointers { virtual void min_xf_pointers(int, double **, double **) {} virtual void min_xf_get(int) {} virtual void min_x_set(int) {} + virtual void transfer_history(double *, double*) {} // management of callbacks to be run from ev_tally() @@ -202,6 +205,7 @@ class Pair : protected Pointers { double tabinner; // inner cutoff for Coulomb table double tabinner_disp; // inner cutoff for dispersion table + public: // custom data type for accessing Coulomb tables