diff --git a/doc/src/pair_granular.txt b/doc/src/pair_granular.txt index 35b64bf29d..8bf4f20d5d 100644 --- a/doc/src/pair_granular.txt +++ b/doc/src/pair_granular.txt @@ -1,3 +1,10 @@ + + + "LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c :link(lws,http://lammps.sandia.gov) @@ -14,323 +21,185 @@ pair_style granular/multi command :h3 pair_style style cutoff :pre style = {granular} or {granular/multi} :ulb,l -cutoff = global cutoff (optional). See discussion below. +cutoff = global cutoff (optional). See discussion below. :l +:ule [Examples:] pair_style granular -pair coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall -pair coeff 2 2 hertz 200.0 20.0 tangential mindlin 300.0 50.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall +pair_coeff * * hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.4 :pre + +pair_style granular +pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall +pair_coeff 2 2 hertz 200.0 20.0 tangential mindlin 300.0 50.0 0.1 rolling sds 200.0 100.0 0.1 twisting marshall :pre pair_style granular/multi -pair coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall -pair coeff 2 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall -pair coeff 1 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall +pair_coeff 1 1 hertz 1000.0 50.0 tangential mindlin 800.0 50.0 0.5 rolling sds 500.0 200.0 0.5 twisting marshall +pair_coeff 2 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall +pair_coeff 1 2 dmt 1000.0 50.0 800.0 10.0 tangential mindlin 800.0 50.0 0.1 roll sds 500.0 200.0 0.1 twisting marshall :pre [Description:] The {granular} styles support a variety of options for the normal, tangential, rolling and twisting -forces resulting from contact between two granular particles. The computed force depends on the combination -of choices for these models. +forces resulting from contact between two granular particles. This expands on the options offered +by the "pair gran/*"_pair_gran.html options. The total computed forces and torques depend on the combination +of choices for these various modes of motion. All model options and parameters are entered in the "pair_coeff"_pair_coeff.html command, as described below. Unlike e.g. "pair gran/hooke"_pair_gran.html, coefficient values are not global, but can be set to different values for various combinations of particle types, as determined by the "pair_coeff"_pair_coeff.html command. -In the case of {granular}, coefficients -can vary between particle types, but model choices cannot. For instance, in the first -example above, the stiffness, damping, and tangential friction are different for type 1 - type 1 and type 2 - type 2 interactions, but -both 1-1 and 2-2 interactions must have the same model form, hence all keywords are identical between the two types. Cross-coefficients -for 1-2 interactions for the case of the {hertz} model above are set via simple geometric mixing rules. The {granular/multi} +For {pair_style granular}, coefficients can vary between particle types, but model choices +cannot. For instance, in the first +example above, the stiffness, damping, and tangential friction are different for +type 1 - type 1 and type 2 - type 2 interactions, but +both 1-1 and 2-2 interactions must have the same model form, hence all keywords are +identical between the two types. Cross-coefficients +for 1-2 interactions for the case of the {hertz} model above are set via simple +geometric mixing rules. The {granular/multi} style removes this restriction at a small cost in computational efficiency, so that different particle types can potentially interact via different model forms. As shown in the second example, 1-1 interactions are based on a Hertzian contact model and 2-2 interactions are based on a {dmt} model (see below). In the case that 1-1 and 2-2 interactions have different model forms, mixing of coefficients cannot be determined, so 1-2 interactions must be explicitly defined via the pair coeff command, otherwise an error results. -The first required keyword for the pair coeff command is the normal contact model. Currently supported options and -the required arguments are: +:line -{hooke} : k_n, damping -{hertz} : k_n, damping +The first required keyword for the {pair_coeff} command is the normal contact model. Currently supported options +for normal contact models and their required arguments are: + +{hooke} : \(k_n\), damping +{hertz} : \(k_n\), damping {hertz/material} : E, damping, G {dmt} : E, damping, G, cohesion -{jkr} : E, damping, G, cohesion +{jkr} : E, damping, G, cohesion :ol -Here, k_n is spring stiffness, damping is a damping constant or a coefficient of restitution, depending on +Here, \(k_n\) is spring stiffness, damping is a damping constant or a coefficient of restitution, depending on the choice of damping model (see below), E and G are Young's modulus and shear modulus, in units of pressure, and cohesion is a surface energy density, in units of energy/length^2. -For the {hooke} model, the normal component of force is given by: -:c,image(Eqs/hooke_normal.jpg) +For the {hooke} model, the normal (elastic) component of force between two particles {i} and {j} is given by: +\begin\{equation\} +\mathbf\{F\}_\{ne, Hooke\} = k_N \delta_\{ij\} \mathbf\{n\} +\end\{equation\} -For {hertz}, the normal force is given by: -:c,image{Eqs/hertz_normal.jpg} +Where \(\delta = R_i + R_j - \|\mathbf\{r\}_\{ij\}\|\) is the particle overlap, +\(R_i, R_j\) are the particle radii, +\(\mathbf\{r\}_\{ij\} = \mathbf\{r\}_j - \mathbf\{r\}_i\) is the vector separating the +two particle centers +and \(\mathbf\{n\} = \frac\{\mathbf\{r\}_\{ij\}\}\{\|\mathbf\{r\}_\{ij\}\|\}\). -For both [hooke] and [hertz], stiffness for unspecified cross-terms is given by simple geometric mixing -(e.g. if stiffness is specified for type 1 and type 2 particles as k_n_1 and k_n_2, respectively, -type 1 - type 2 contacts use a stiffness given by k_n_{12} = sqrt(k_n_1*k_n_2)) +For the {hertz} model, the normal component of force is given by: +\begin\{equation\} +\mathbf\{F\}_\{ne, Hertz\} = k_N R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\} +\end\{equation\} -For {hertz/material}, the form is the same as above, but coefficients are computed differently, and mixing follows -a different rule based on shear modulus: -:c,image{Eqs/hertz_material_normal.jpg} +Here, \(R_\{eff\} = \frac\{R_i R_j\}\{R_i + R_j\}\) is the effective radius, denoted for simplicity as {R} from here on. -For {dmt}, the normal force is given by: -:c,image{Eqs/dmt_normal.jpg} +For the {hertz/material} model, the force is given by: +\begin\{equation\} +\mathbf\{F\}_\{ne, Hertz/material\} = \frac\{4\}\{3\} E_\{eff\} R_\{eff\}^\{1/2\}\delta_\{ij\}^\{3/2\} \mathbf\{n\} +\end\{equation\} -Where gamma is cohesion. +Here, \(E_\{eff\} = E = \left(\frac\{1-\nu_i^2\}\{E_i\} + \frac\{1-\nu_j^2\}\{E_j\}\right)^\{-1\}\) +is the effectve Young's modulus, +with \(\nu_i, \nu_j \) the Poisson ratios of the particles, which are related to the +input shear and Young's moduli by \(\nu_i = E_i/2G_i - 1\). Thus, if the elastic and shear moduli of the +two particles are the same, the {hertz/material} +model is equivalent to the {hertz} model with \(k_N = 4/3 E_\{eff\}\) -For {jkr}, the normal force is given by: -:c,image{Eqs/jkr_normal.jpg} +The {dmt} model corresponds to the Derjaguin-Muller-Toporov model, +where the force is simply Hertz with an additional attractive cohesion term: +\begin\{equation\} +\mathbf\{F\}_\{ne, dmt\} = \left(\frac\{4\}\{3\} E R^\{1/2\}\delta_\{ij\}^\{3/2\} - 4\pi\gamma R\right)\mathbf\{n\} +\end\{equation\} -The same mixing rule for stiffness as for {hertz/material} is used by both the {dmt} and {jkr} models. +The {jkr} model is the Johnson-Kendall-Roberts model, where the force is computed as: +\begin\{equation\} +\label\{eq:force_jkr\} +\mathbf\{F\}_\{ne, jkr\} = \left(\frac\{4Ea^3\}\{3R\} - 2\pi a^2\sqrt\{\frac\{4\gamma E\}\{\pi a\}\}\right)\mathbf\{n\} +\end\{equation\} + +Here, {a} is the radius of the contact zone, related to the overlap \(\delta\) according to: +\begin\{equation\} +\delta = a^2/R - 2\sqrt\{\pi \gamma a/E\} +\end\{equation\} -The tangential contact model must also be specified, which follows -the required {tangential} keyword. Currently supported options -and their required arguments are: +LAMMPS internally inverts the equation above to solve for {a} in terms of \(\delta\), then solves for +the force in the previous equation. Additionally, note that the JKR model allows for a tensile force beyond +contact (i.e. for \(\delta < 0\)), up to a maximum tensile force of \(-3\pi\gamma R\) (also known as +the 'pull-off' force). +Note that this is a hysteretic effect, where particles that are not contacting initially +will not experience force until they come into contact \(\delta \geq 0\); as they move apart +and (\(\delta < 0\)), they experience a tensile force up to \(-3\pi\gamma R\), +at which point they will lose contact. -{no_history}: k_t, tangential_damping, friction coefficient -{mindlin}: k_t, tangential_damping, friction coefficient +In addition to the above options, the normal force is augmented by a damping term. The optional +{damping} keyword to the {pair_coeff} command followed by the model choice determines the form of the damping. +The damping coefficient that was specified for the normal model +settings is used in computing the damping term, as described below. Note this damping parameter +may be interpreted differently depending on the model choice. +The options for the damping model currently supported are: + +{velocity} +{viscoelastic} +{tsuji} :ol + +If the {damping} keyword is not specified, the {viscoelastic} model is used by default. + +For {damping velocity}, the normal damping is simply proportional to the velocity: +\begin\{equation\} +F_\{N,damp\} = -\gamma_N\mathbf\{v\}_\{N,rel\} +\end\{equation\} + +Here, \(\gamma_N\) is the damping coefficient, in units of {mass}/{time}, +\(\mathbf\{v\}_\{N,rel\} = (\mathbf\{v\}_i - \mathbf\{v\}_j) \cdot \mathbf\{n\}\) +is the component of relative velocity along the direction of the vector \(\mathbf\{n\}\) that connects the centers of +particles {i} and {j}. + +The {damping viscoelastic} model is based on the viscoelastic treatment of "(Brilliantov et al)"_#Brill1996, +where the normal damping is given by: +\begin\{equation\} +F_\{N,damp\} = -\gamma_N a m_\{eff\} \mathbf\{v\}_\{N,rel\} +\end\{equation\} + +Here, \(m_\{eff\} = m_i m_j/(m_i + m_j)\) is the effective mass, {a} is the contact radius, given by \(a =\sqrt\{R\delta\}\) +for all models except {jkr}, for which it is given implicitly according to \(delta = a^2/R - 2\sqrt\{\pi \gamma a/E\}\). +In this case, \(\gamma_N\) is the damping coefficient, in units of 1/({time}*{distance}). + +The {tsuji} model is based on the work of "(Tsuji et al)"_#Tsuji1992. Here, the + + :line + +Following the normal contact model settings, the {pair_coeff} command requires specification +of the tangential contact model. The required keyword {tangential} is expected, followed by the model choice and associated +parameters. Currently supported tangential model choices and their expected parameters are as follows: -For {no_history}, the tangential force is computed according to: -:c,image{Eqs/tangential_nohistory.jpg} +{nohistory} : \(\gamma_t\), \(\mu_s\) +{history} : \(k_t\), \(\gamma_t\), \(\mu_s\) :ol -For {mindlin}, tangential force is: -:c,image{Eqs/tangential_mindlin.jpg} +Here, \(\gamma_t\) is the tangential damping coefficient, \(\mu_s\) is the tangential (or sliding) friction +coefficient, and \(k_t\) is the tangential stiffness. -The total force on a particle is the sum of the normal and tangential forces from all interactions. The tangential -force also induces a torque on both particles in a contacting pair. Additionally, rolling and twisting friction -models can also be applied, which may induce additional torques (but no force). The following options are -supported for the rolling friction model +For {nohistory}, a simple velocity-dependent Coulomb friction criterion is used, which reproduces the behavior +of the {pair gran/hooke} style. The tangential force (\mathbf\{F\}_t\) is given by: +\begin\{equation\} +\mathbf\{F\}_t = -min(\mu_s \|\mathbf\{F\}_n\|, \gamma_t m_\{eff\}\|\mathbf\{v\}_\{t, rel\}\|) \mathbf\{t\} +\end\{equation\} + +Where \(\|\mathbf\{F\}_n\) is the magnitude of the normal force, +\(\mathbf\{v\}_\{t, rel\} = \mathbf\{v\}_\{t\} - (R_i\Omega_i + R_j\Omega_j) \times \mathbf\{n\}\) is the relative tangential +velocity at the point of contact, \(\mathbf\{v\}_\{t\} = \mathbf\{v\}_n - \) + + + :link(Brill1996) +[(Brilliantov et al, 1996)] Brilliantov, N. V., Spahn, F., Hertzsch, J. M., & Poschel, T. (1996). +Model for collisions in granular gases. Physical review E, 53(5), 5382. -The first required keyword -in the pair coeff command is the choice -of normal force contact model, for which current opitons are {hooke}, {hertz} + :link(Tsuji1992) + [(Tsuji et al, 1992)] Tsuji, Y., Tanaka, T., & Ishida, T. (1992). Lagrangian numerical simulation of plug flow of + cohesionless particles in a horizontal pipe. Powder technology, 71(3), 239-250. -The {gran} styles use the following formulas for the frictional force -between two granular particles, as described in -"(Brilliantov)"_#Brilliantov, "(Silbert)"_#Silbert, and -"(Zhang)"_#Zhang3, when the distance r between two particles of radii -Ri and Rj is less than their contact distance d = Ri + Rj. There is -no force between the particles when r > d. - -The two Hookean styles use this formula: - -:c,image(Eqs/pair_gran_hooke.jpg) - -The Hertzian style uses this formula: - -:c,image(Eqs/pair_gran_hertz.jpg) - -In both equations the first parenthesized term is the normal force -between the two particles and the second parenthesized term is the -tangential force. The normal force has 2 terms, a contact force and a -damping force. The tangential force also has 2 terms: a shear force -and a damping force. The shear force is a "history" effect that -accounts for the tangential displacement between the particles for the -duration of the time they are in contact. This term is included in -pair styles {hooke/history} and {hertz/history}, but is not included -in pair style {hooke}. The tangential damping force term is included -in all three pair styles if {dampflag} is set to 1; it is not included -if {dampflag} is set to 0. - -The other quantities in the equations are as follows: - -delta = d - r = overlap distance of 2 particles -Kn = elastic constant for normal contact -Kt = elastic constant for tangential contact -gamma_n = viscoelastic damping constant for normal contact -gamma_t = viscoelastic damping constant for tangential contact -m_eff = Mi Mj / (Mi + Mj) = effective mass of 2 particles of mass Mi and Mj -Delta St = tangential displacement vector between 2 particles \ - which is truncated to satisfy a frictional yield criterion -n_ij = unit vector along the line connecting the centers of the 2 particles -Vn = normal component of the relative velocity of the 2 particles -Vt = tangential component of the relative velocity of the 2 particles :ul - -The Kn, Kt, gamma_n, and gamma_t coefficients are specified as -parameters to the pair_style command. If a NULL is used for Kt, then -a default value is used where Kt = 2/7 Kn. If a NULL is used for -gamma_t, then a default value is used where gamma_t = 1/2 gamma_n. - -The interpretation and units for these 4 coefficients are different in -the Hookean versus Hertzian equations. - -The Hookean model is one where the normal push-back force for two -overlapping particles is a linear function of the overlap distance. -Thus the specified Kn is in units of (force/distance). Note that this -push-back force is independent of absolute particle size (in the -monodisperse case) and of the relative sizes of the two particles (in -the polydisperse case). This model also applies to the other terms in -the force equation so that the specified gamma_n is in units of -(1/time), Kt is in units of (force/distance), and gamma_t is in units -of (1/time). - -The Hertzian model is one where the normal push-back force for two -overlapping particles is proportional to the area of overlap of the -two particles, and is thus a non-linear function of overlap distance. -Thus Kn has units of force per area and is thus specified in units of -(pressure). The effects of absolute particle size (monodispersity) -and relative size (polydispersity) are captured in the radii-dependent -pre-factors. When these pre-factors are carried through to the other -terms in the force equation it means that the specified gamma_n is in -units of (1/(time*distance)), Kt is in units of (pressure), and -gamma_t is in units of (1/(time*distance)). - -Note that in the Hookean case, Kn can be thought of as a linear spring -constant with units of force/distance. In the Hertzian case, Kn is -like a non-linear spring constant with units of force/area or -pressure, and as shown in the "(Zhang)"_#Zhang3 paper, Kn = 4G / -(3(1-nu)) where nu = the Poisson ratio, G = shear modulus = E / -(2(1+nu)), and E = Young's modulus. Similarly, Kt = 4G / (2-nu). -(NOTE: in an earlier version of the manual, we incorrectly stated that -Kt = 8G / (2-nu).) - -Thus in the Hertzian case Kn and Kt can be set to values that -corresponds to properties of the material being modeled. This is also -true in the Hookean case, except that a spring constant must be chosen -that is appropriate for the absolute size of particles in the model. -Since relative particle sizes are not accounted for, the Hookean -styles may not be a suitable model for polydisperse systems. - -NOTE: In versions of LAMMPS before 9Jan09, the equation for Hertzian -interactions did not include the sqrt(RiRj/Ri+Rj) term and thus was -not as accurate for polydisperse systems. For monodisperse systems, -sqrt(RiRj/Ri+Rj) is a constant factor that effectively scales all 4 -coefficients: Kn, Kt, gamma_n, gamma_t. Thus you can set the values -of these 4 coefficients appropriately in the current code to reproduce -the results of a previous Hertzian monodisperse calculation. For -example, for the common case of a monodisperse system with particles -of diameter 1, all 4 of these coefficients should now be set 2x larger -than they were previously. - -Xmu is also specified in the pair_style command and is the upper limit -of the tangential force through the Coulomb criterion Ft = xmu*Fn, -where Ft and Fn are the total tangential and normal force components -in the formulas above. Thus in the Hookean case, the tangential force -between 2 particles grows according to a tangential spring and -dash-pot model until Ft/Fn = xmu and is then held at Ft = Fn*xmu until -the particles lose contact. In the Hertzian case, a similar analogy -holds, though the spring is no longer linear. - -NOTE: Normally, xmu should be specified as a fractional value between -0.0 and 1.0, however LAMMPS allows large values (up to 1.0e4) to allow -for modeling of systems which can sustain very large tangential -forces. - -The effective mass {m_eff} is given by the formula above for two -isolated particles. If either particle is part of a rigid body, its -mass is replaced by the mass of the rigid body in the formula above. -This is determined by searching for a "fix rigid"_fix_rigid.html -command (or its variants). - -For granular styles there are no additional coefficients to set for -each pair of atom types via the "pair_coeff"_pair_coeff.html command. -All settings are global and are made via the pair_style command. -However you must still use the "pair_coeff"_pair_coeff.html for all -pairs of granular atom types. For example the command - -pair_coeff * * :pre - -should be used if all atoms in the simulation interact via a granular -potential (i.e. one of the pair styles above is used). If a granular -potential is used as a sub-style of "pair_style -hybrid"_pair_hybrid.html, then specific atom types can be used in the -pair_coeff command to determine which atoms interact via a granular -potential. - -:line - -Styles with a {gpu}, {intel}, {kk}, {omp}, or {opt} suffix are -functionally the same as the corresponding style without the suffix. -They have been optimized to run faster, depending on your available -hardware, as discussed on the "Speed packages"_Speed_packages.html doc -page. The accelerated styles take the same arguments and should -produce the same results, except for round-off and precision issues. - -These accelerated styles are part of the GPU, USER-INTEL, KOKKOS, -USER-OMP and OPT packages, respectively. They are only enabled if -LAMMPS was built with those packages. See the "Build -package"_Build_package.html doc page for more info. - -You can specify the accelerated styles explicitly in your input script -by including their suffix, or you can use the "-suffix command-line -switch"_Run_options.html when you invoke LAMMPS, or you can use the -"suffix"_suffix.html command in your input script. - -See the "Speed packages"_Speed_packages.html doc page for more -instructions on how to use the accelerated styles effectively. - -:line - -[Mixing, shift, table, tail correction, restart, rRESPA info]: - -The "pair_modify"_pair_modify.html mix, shift, table, and tail options -are not relevant for granular pair styles. - -These pair styles write their information to "binary restart -files"_restart.html, so a pair_style command does not need to be -specified in an input script that reads a restart file. - -These pair styles can only be used via the {pair} keyword of the -"run_style respa"_run_style.html command. They do not support the -{inner}, {middle}, {outer} keywords. - -The single() function of these pair styles returns 0.0 for the energy -of a pairwise interaction, since energy is not conserved in these -dissipative potentials. It also returns only the normal component of -the pairwise interaction force. However, the single() function also -calculates 10 extra pairwise quantities. The first 3 are the -components of the tangential force between particles I and J, acting -on particle I. The 4th is the magnitude of this tangential force. -The next 3 (5-7) are the components of the relative velocity in the -normal direction (along the line joining the 2 sphere centers). The -last 3 (8-10) the components of the relative velocity in the -tangential direction. - -These extra quantities can be accessed by the "compute -pair/local"_compute_pair_local.html command, as {p1}, {p2}, ..., -{p10}. - -:line - -[Restrictions:] - -All the granular pair styles are part of the GRANULAR package. It is -only enabled if LAMMPS was built with that package. See the "Build -package"_Build_package.html doc page for more info. - -These pair styles require that atoms store torque and angular velocity -(omega) as defined by the "atom_style"_atom_style.html. They also -require a per-particle radius is stored. The {sphere} atom style does -all of this. - -This pair style requires you to use the "comm_modify vel -yes"_comm_modify.html command so that velocities are stored by ghost -atoms. - -These pair styles will not restart exactly when using the -"read_restart"_read_restart.html command, though they should provide -statistically similar results. This is because the forces they -compute depend on atom velocities. See the -"read_restart"_read_restart.html command for more details. - -[Related commands:] - -"pair_coeff"_pair_coeff.html - -[Default:] none - -:line - -:link(Brilliantov) -[(Brilliantov)] Brilliantov, Spahn, Hertzsch, Poschel, Phys Rev E, 53, -p 5382-5392 (1996). - -:link(Silbert) -[(Silbert)] Silbert, Ertas, Grest, Halsey, Levine, Plimpton, Phys Rev -E, 64, p 051302 (2001). - -:link(Zhang3) -[(Zhang)] Zhang and Makse, Phys Rev E, 72, p 011301 (2005). + + \ No newline at end of file diff --git a/doc/src/pairs.txt b/doc/src/pairs.txt index ca79051053..3cd20c728d 100644 --- a/doc/src/pairs.txt +++ b/doc/src/pairs.txt @@ -41,6 +41,7 @@ Pair Styles :h1 pair_gauss pair_gayberne pair_gran + pair_granular pair_gromacs pair_gw pair_hbond_dreiding diff --git a/src/GRANULAR/pair_granular.cpp b/src/GRANULAR/pair_granular.cpp index 82a470f83b..3713b9251c 100644 --- a/src/GRANULAR/pair_granular.cpp +++ b/src/GRANULAR/pair_granular.cpp @@ -11,8 +11,9 @@ See the README file in the top-level LAMMPS directory. ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- -Contributing authors: Leo Silbert (SNL), Gary Grest (SNL), - Jeremy Lechman (SNL), Dan Bolintineanu (SNL), Ishan Srivastava (SNL) +Contributing authors: +Dan Bolintineanu (SNL), Ishan Srivastava (SNL), Jeremy Lechman(SNL) +Leo Silbert (SNL), Gary Grest (SNL) ----------------------------------------------------------------------- */ #include @@ -50,7 +51,7 @@ using namespace MathConst; enum {HOOKE, HERTZ, HERTZ_MATERIAL, DMT, JKR}; enum {VELOCITY, VISCOELASTIC, TSUJI}; -enum {TANGENTIAL_NOHISTORY, TANGENTIAL_MINDLIN}; +enum {TANGENTIAL_NOHISTORY, TANGENTIAL_HISTORY, TANGENTIAL_MINDLIN}; enum {TWIST_NONE, TWIST_NOHISTORY, TWIST_SDS, TWIST_MARSHALL}; enum {ROLL_NONE, ROLL_NOHISTORY, ROLL_SDS}; @@ -98,13 +99,24 @@ PairGranular::~PairGranular() if (allocated) { memory->destroy(setflag); memory->destroy(cutsq); - memory->destroy(cut); + memory->destroy(cutoff_type); memory->destroy(normal_coeffs); memory->destroy(tangential_coeffs); memory->destroy(roll_coeffs); memory->destroy(twist_coeffs); + memory->destroy(Emod); + memory->destroy(poiss); + + memory->destroy(normal_model); + memory->destroy(damping_model); + memory->destroy(tangential_model); + memory->destroy(roll_model); + memory->destroy(twist_model); + + + delete [] onerad_dynamic; delete [] onerad_frozen; delete [] maxrad_dynamic; @@ -113,725 +125,7 @@ PairGranular::~PairGranular() memory->destroy(mass_rigid); } -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_twist, int Tp_roll > -void PairGranular::compute_templated(int eflag, int vflag) -#else -void PairGranular::compute_untemplated - (int Tp_normal, int Tp_damping, int Tp_tangential, - int Tp_twist, int Tp_roll, int eflag, int vflag) -#endif +void PairGranular::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; @@ -847,14 +141,14 @@ void PairGranular::compute_untemplated double Fne, Ft, Fdamp, Fntot, Fcrit, Fscrit, Frcrit; double fs, fs1, fs2, fs3; + double mi,mj,meff,damp,ccel,tor1,tor2,tor3; + double relrot1,relrot2,relrot3,vrl1,vrl2,vrl3,vrlmag,vrlmaginv; + //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; @@ -947,7 +241,7 @@ void PairGranular::compute_untemplated Reff = radi*radj/(radi+radj); touchflag = false; - if (Tp_normal == JKR){ + if (normal_model[itype][jtype] == JKR){ if (touch[jj]){ R2 = Reff*Reff; coh = normal_coeffs[itype][jtype][3]; @@ -1008,7 +302,7 @@ void PairGranular::compute_untemplated delta = radsum - r; dR = delta*Reff; - if (Tp_normal == JKR){ + if (normal_model[itype][jtype] == JKR){ touch[jj] = 1; R2=Reff*Reff; coh = normal_coeffs[itype][jtype][3]; @@ -1031,22 +325,21 @@ void PairGranular::compute_untemplated else{ knfac = E; //Hooke Fne = knfac*delta; - if (Tp_normal != HOOKE) - a = sqrt(dR); + a = sqrt(dR); + if (normal_model[itype][jtype] != HOOKE) Fne *= a; - if (Tp_normal == DMT) + if (normal_model[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 (Tp_damping == VELOCITY){ + if (damping_model[itype][jtype] == VELOCITY){ damp_normal = 1; } - else if (Tp_damping == VISCOELASTIC){ - if (Tp_normal == HOOKE) a = sqrt(dR); + else if (damping_model[itype][jtype] == VISCOELASTIC){ damp_normal = a*meff; } - else if (Tp_damping == TSUJI){ + else if (damping_model[itype][jtype] == TSUJI){ damp_normal = sqrt(meff*knfac); } @@ -1081,8 +374,7 @@ void PairGranular::compute_untemplated history = &allhistory[size_history*jj]; } - - if (Tp_normal == JKR){ + if (normal_model[itype][jtype] == JKR){ F_pulloff = 3*M_PI*coh*Reff; Fcrit = fabs(Fne + 2*F_pulloff); } @@ -1096,7 +388,7 @@ void PairGranular::compute_untemplated k_tangential = tangential_coeffs[itype][jtype][0]; damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; - if (Tp_tangential > 0){ + if (tangential_history){ shrmag = sqrt(history[0]*history[0] + history[1]*history[1] + history[2]*history[2]); @@ -1153,7 +445,7 @@ void PairGranular::compute_untemplated // Rolling resistance //**************************************** - if (Tp_roll != ROLL_NONE){ + if (roll_model[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]; @@ -1168,7 +460,7 @@ void PairGranular::compute_untemplated if (vrlmag != 0.0) vrlmaginv = 1.0/vrlmag; else vrlmaginv = 0.0; - if (Tp_roll > 1){ + if (roll_history){ int rhist0 = roll_history_index; int rhist1 = rhist0 + 1; int rhist2 = rhist1 + 1; @@ -1197,6 +489,7 @@ void PairGranular::compute_untemplated 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; @@ -1231,9 +524,9 @@ void PairGranular::compute_untemplated //**************************************** // Twisting torque, including history effects //**************************************** - if (Tp_twist != TWIST_NONE){ + if (twist_model[itype][jtype] != TWIST_NONE){ magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (Tp_twist == TWIST_MARSHALL){ + if (twist_model[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; @@ -1243,7 +536,7 @@ void PairGranular::compute_untemplated damp_twist = twist_coeffs[itype][jtype][1]; mu_twist = twist_coeffs[itype][jtype][2]; } - if (Tp_twist > 1){ + if (twist_model[itype][jtype] > 1){ if (historyupdate){ history[twist_history_index] += magtwist*dt; } @@ -1278,7 +571,7 @@ void PairGranular::compute_untemplated torque[i][1] -= radi*tor2; torque[i][2] -= radi*tor3; - if (Tp_twist != TWIST_NONE){ + if (twist_model[itype][jtype] != TWIST_NONE){ tortwist1 = magtortwist * nx; tortwist2 = magtortwist * ny; tortwist3 = magtortwist * nz; @@ -1288,7 +581,7 @@ void PairGranular::compute_untemplated torque[i][2] += tortwist3; } - if (Tp_roll != ROLL_NONE){ + if (roll_model[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); @@ -1307,12 +600,12 @@ void PairGranular::compute_untemplated torque[j][1] -= radj*tor2; torque[j][2] -= radj*tor3; - if (Tp_twist != TWIST_NONE){ + if (twist_model[itype][jtype] != TWIST_NONE){ torque[j][0] -= tortwist1; torque[j][1] -= tortwist2; torque[j][2] -= tortwist3; } - if (Tp_roll != ROLL_NONE){ + if (roll_model[itype][jtype] != ROLL_NONE){ torque[j][0] -= torroll1; torque[j][1] -= torroll2; torque[j][2] -= torroll3; @@ -1341,12 +634,21 @@ void PairGranular::allocate() setflag[i][j] = 0; memory->create(cutsq,n+1,n+1,"pair:cutsq"); - memory->create(cut,n+1,n+1,"pair:cut"); + memory->create(cutoff_type,n+1,n+1,"pair:cutoff_type"); 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(Emod,n+1,n+1,"pair:Emod"); + memory->create(poiss,n+1,n+1,"pair:poiss"); + + memory->create(normal_model,n+1,n+1,"pair:normal_model"); + memory->create(damping_model,n+1,n+1,"pair:damping_model"); + memory->create(tangential_model,n+1,n+1,"pair:tangential_model"); + memory->create(roll_model,n+1,n+1,"pair:roll_model"); + memory->create(twist_model,n+1,n+1,"pair:twist_model"); + onerad_dynamic = new double[n+1]; onerad_frozen = new double[n+1]; maxrad_dynamic = new double[n+1]; @@ -1365,9 +667,9 @@ void PairGranular::settings(int narg, char **arg) else{ cutoff_global = -1; //Will be set based on particle sizes, model choice } - tangential_history = 0; + + normal_history = tangential_history = 0; roll_history = twist_history = 0; - normal_set = tangential_set = damping_set = roll_set = twist_set = 0; } /* ---------------------------------------------------------------------- @@ -1376,10 +678,15 @@ void PairGranular::settings(int narg, char **arg) void PairGranular::coeff(int narg, char **arg) { - double normal_coeffs_local[4]; - double tangential_coeffs_local[4]; - double roll_coeffs_local[4]; - double twist_coeffs_local[4]; + int normal_model_one, damping_model_one; + int tangential_model_one, roll_model_one, twist_model_one; + + double normal_coeffs_one[4]; + double tangential_coeffs_one[4]; + double roll_coeffs_one[4]; + double twist_coeffs_one[4]; + + double cutoff_one = -1; if (narg < 2) error->all(FLERR,"Incorrect args for pair coefficients"); @@ -1390,209 +697,200 @@ 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); + //Defaults + normal_model_one = tangential_model_one = -1; + roll_model_one = twist_model_one = 0; + damping_model_one = VISCOELASTIC; + 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"); - if (!normal_set) normal = HOOKE; - else if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); - 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; + normal_model_one = HOOKE; + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs_one[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"); - if (!normal_set) normal = HERTZ; - else if (normal_set && normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); - 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; + normal_model_one = HERTZ; + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //kn + normal_coeffs_one[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"); - if (!normal_set) normal = HERTZ_MATERIAL; - else if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); - 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; + normal_model_one = HERTZ; + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio 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"); - if (!normal_set) normal = DMT; - else if (normal != DMT) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); - 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; + normal_model_one = DMT; + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1])*FOURTHIRDS; //E + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio + normal_coeffs_one[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; - if (!normal_set) normal = JKR; - else if (normal != JKR) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); - 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; + normal_model_one = JKR; + normal_coeffs_one[0] = force->numeric(FLERR,arg[iarg+1]); //E + normal_coeffs_one[1] = force->numeric(FLERR,arg[iarg+2]); //damping + normal_coeffs_one[2] = force->numeric(FLERR,arg[iarg+3]); //Poisson's ratio + normal_coeffs_one[3] = force->numeric(FLERR,arg[iarg+4]); //cohesion iarg += 5; } else if (strcmp(arg[iarg], "damp") == 0){ if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model"); if (strcmp(arg[iarg+1], "velocity") == 0){ - if (!damping_set) damping = VELOCITY; - else if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + damping_model_one = VELOCITY; + iarg += 1; } else if (strcmp(arg[iarg+1], "viscoelastic") == 0){ - if (!damping_set) damping = VISCOELASTIC; - else if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + damping_model_one = VISCOELASTIC; + iarg += 1; } - else if (strcmp(arg[iarg+1], "tsuji") == 0){ - if (!damping_set) damping = TSUJI; - if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + else if (strcmp(arg[iarg], "tsuji") == 0){ + damping_model_one = TSUJI; + iarg += 1; } - damping_set = 1; - iarg += 2; } else if (strcmp(arg[iarg], "tangential") == 0){ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); - if (strcmp(arg[iarg+1], "nohistory") == 0){ - if (!tangential_set) tangential = TANGENTIAL_NOHISTORY; - else if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types"); + if (strcmp(arg[iarg+1], "linear_nohistory") == 0){ + tangential_model_one = TANGENTIAL_NOHISTORY; } - else if (strcmp(arg[iarg+1], "mindlin") == 0){ - if (!tangential_set) tangential = TANGENTIAL_MINDLIN; - else if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types");; + else if (strcmp(arg[iarg+1], "linear_history") == 0){ + tangential_model_one = TANGENTIAL_HISTORY; 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. + tangential_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt + tangential_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + tangential_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. iarg += 5; } else if (strcmp(arg[iarg], "rolling") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); if (strcmp(arg[iarg+1], "none") == 0){ - if (!roll_set) roll = ROLL_NONE; - else if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); + roll_model_one = ROLL_NONE; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); if (strcmp(arg[iarg+1], "nohistory") == 0){ - if (!roll_set) roll = ROLL_NOHISTORY; - else if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); + roll_model_one = ROLL_NOHISTORY; } else if (strcmp(arg[iarg+1], "sds") == 0){ - if (!roll_set) roll = ROLL_SDS; - else if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); + roll_model_one = ROLL_SDS; 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]); //kR - roll_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR - roll_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. + roll_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kR + roll_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR + roll_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. iarg += 5; } } else if (strcmp(arg[iarg], "twisting") == 0){ if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); if (strcmp(arg[iarg+1], "none") == 0){ - if (!twist_set) twist = TWIST_NONE; - else if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + twist_model_one = TWIST_NONE; iarg += 2; } else if (strcmp(arg[iarg+1], "marshall") == 0){ - if (!twist_set) twist = TWIST_MARSHALL; - else if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + twist_model_one = TWIST_MARSHALL; twist_history = 1; - twist_set = 1; iarg += 2; } else{ if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); - else if (strcmp(arg[iarg+1], "nohistory") == 0){ - if (!twist_set) twist = TWIST_NOHISTORY; - if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + if (strcmp(arg[iarg+1], "nohistory") == 0){ + twist_model_one = TWIST_NOHISTORY; } else if (strcmp(arg[iarg+1], "sds") == 0){ - if (!twist_set) twist = TWIST_SDS; - else if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + twist_model_one = TWIST_SDS; 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. + twist_coeffs_one[0] = force->numeric(FLERR,arg[iarg+2]); //kt + twist_coeffs_one[1] = force->numeric(FLERR,arg[iarg+3]); //gammat + twist_coeffs_one[2] = force->numeric(FLERR,arg[iarg+4]); //friction coeff. iarg += 5; } } + else if (strcmp(arg[iarg], "cutoff") == 0){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + cutoff_one = force->numeric(FLERR,arg[iarg+1]); + } else error->all(FLERR, "Illegal pair coeff command"); } //It is an error not to specify normal or tangential model - if (!normal_set) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model"); - if (!tangential_set) error->all(FLERR, "Illegal pair coeff command, must specify tangential contact model"); - - //If unspecified, set damping to VISCOELASTIC, twist/roll to NONE (cannot be changed by subsequent pair_coeff commands) - if (!damping_set) damping = VISCOELASTIC; - if (!roll_set) roll = ROLL_NONE; - if (!twist_set) twist = TWIST_NONE; - damping_set = roll_set = twist_set = 1; + if ((normal_model_one < 0) || (tangential_model_one < 0)) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model"); int count = 0; double damp; - if (damping == TSUJI){ + if (damping_model_one == TSUJI){ double cor; - cor = normal_coeffs_local[1]; + cor = normal_coeffs_one[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]; + else damp = normal_coeffs_one[1]; for (int i = ilo; i <= ihi; i++) { for (int j = MAX(jlo,i); j <= jhi; j++) { - normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_local[0]; + normal_model[i][j] = normal_model[j][i] = normal_model_one; normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp; - 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[j][i][3] = normal_coeffs_local[3]; + if (normal_model_one != HERTZ && normal_model_one != HOOKE){ + Emod[i][j] = Emod[j][i] = normal_coeffs_one[0]; + poiss[i][j] = poiss[j][i] = normal_coeffs_one[2]; + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(Emod[i][j], Emod[i][j], poiss[i][j], poiss[i][j]); + } + else{ + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_one[0]; + } + if ((normal_model_one == JKR) || (normal_model_one == DMT)) + normal_coeffs[i][j][3] = normal_coeffs[j][i][3] = normal_coeffs_one[3]; + damping_model[i][j] = damping_model[j][i] = damping_model_one; + + tangential_model[i][j] = tangential_model[j][i] = tangential_model_one; for (int k = 0; k < 3; k++) - tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_local[k]; + tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_one[k]; - if (roll != ROLL_NONE) + roll_model[i][j] = roll_model[j][i] = roll_model_one; + if (roll_model_one != ROLL_NONE) for (int k = 0; k < 3; k++) - roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_local[k]; + roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_one[k]; - if (twist != TWIST_NONE && twist != TWIST_MARSHALL) + twist_model[i][j] = twist_model[j][i] = twist_model_one; + if (twist_model_one != TWIST_NONE && twist_model_one != TWIST_MARSHALL) for (int k = 0; k < 3; k++) - twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_local[k]; + twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_one[k]; + + + cutoff_type[i][j] = cutoff_type[j][i] = cutoff_one; setflag[i][j] = 1; count++; } } - if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); } @@ -1615,7 +913,9 @@ 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 - if (normal == JKR) use_history = 1; + for (int i = 1; i <= atom->ntypes; i++) + for (int j = 1; j <= atom->ntypes; j++) + if (normal_model[i][j] == JKR) use_history = 1; size_history = 3*tangential_history + 3*roll_history + twist_history; @@ -1693,13 +993,11 @@ void PairGranular::init_style() if (ipour >= 0) { itype = i; double radmax = *((double *) modify->fix[ipour]->extract("radius",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 == JKR) radmax = radmax - 0.5*pulloff_distance(radmax, itype); onerad_dynamic[i] = radmax; } } @@ -1711,9 +1009,6 @@ void PairGranular::init_style() for (i = 0; i < nlocal; i++){ double radius_cut = radius[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); } @@ -1743,32 +1038,38 @@ void PairGranular::init_style() double PairGranular::init_one(int i, int j) { double cutoff; - if (setflag[i][j] == 0) { - if (normal != HOOKE && normal != HERTZ){ - normal_coeffs[i][j][0] = normal_coeffs[j][i][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] = normal_coeffs[j][i][2] = mix_stiffnessG(normal_coeffs[i][i][0], normal_coeffs[j][j][0], - normal_coeffs[i][i][2], normal_coeffs[j][j][2]); + if (setflag[i][j] == 0){ + if ((normal_model[i][i] != normal_model[j][j]) || + (damping_model[i][i] != damping_model[j][j]) || + (tangential_model[i][i] != tangential_model[j][j]) || + (roll_model[i][i] != roll_model[j][j]) || + (twist_model[i][i] != twist_model[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); } - else{ + + if (normal_model[i][j] == HERTZ || normal_model[i][j] == HOOKE) normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); - } + else + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = mix_stiffnessE(Emod[i][i], Emod[j][j], poiss[i][i], poiss[j][j]); normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); - if ((normal == JKR) || (normal == DMT)) + if ((normal_model[i][j] == JKR) || (normal_model[i][j] == DMT)) normal_coeffs[i][j][3] = normal_coeffs[j][i][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] = tangential_coeffs[j][i][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); - if (roll != ROLL_NONE){ + if (roll_model[i][j] != ROLL_NONE){ for (int k = 0; k < 3; k++) roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); } - if (twist != TWIST_NONE && twist != TWIST_MARSHALL){ + if (twist_model[i][j] != TWIST_NONE && twist_model[i][j] != TWIST_MARSHALL){ for (int k = 0; k < 3; k++) twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = mix_geom(twist_coeffs[i][i][k], twist_coeffs[j][j][k]); } @@ -1780,22 +1081,40 @@ double PairGranular::init_one(int i, int j) // 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. - + double pulloff; if (cutoff_global < 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]); + if (cutoff_type[i][j] < 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]; + if (normal_model[i][j] == JKR){ + pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_dynamic[j], i, j); + cutoff += pulloff; + } + else{ + pulloff = 0; + } + + if (normal_model[i][j] == JKR) + pulloff = pulloff_distance(maxrad_frozen[i], maxrad_dynamic[j], i, j); + cutoff = MAX(cutoff, maxrad_frozen[i]+maxrad_dynamic[j]+pulloff); + + if (normal_model[i][j] == JKR) + pulloff = pulloff_distance(maxrad_dynamic[i], maxrad_frozen[j], i, j); + cutoff = MAX(cutoff,maxrad_dynamic[i]+maxrad_frozen[j]+pulloff); } - cutoff = cutmax; + 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; + } + } + else{ + cutoff = cutoff_type[i][j]; } } else{ @@ -1804,7 +1123,6 @@ double PairGranular::init_one(int i, int j) return cutoff; } - /* ---------------------------------------------------------------------- proc 0 writes to restart file ------------------------------------------------------------------------- */ @@ -1812,15 +1130,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_model[i][j],sizeof(int),1,fp); + fwrite(&damping_model[i][j],sizeof(int),1,fp); + fwrite(&tangential_model[i][j],sizeof(int),1,fp); + fwrite(&roll_model[i][j],sizeof(int),1,fp); + fwrite(&twist_model[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); @@ -1840,30 +1158,28 @@ 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_model[i][j],sizeof(int),1,fp); + fread(&damping_model[i][j],sizeof(int),1,fp); + fread(&tangential_model[i][j],sizeof(int),1,fp); + fread(&roll_model[i][j],sizeof(int),1,fp); + fread(&twist_model[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_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&damping_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&tangential_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&roll_model[i][j],1,MPI_INT,0,world); + MPI_Bcast(&twist_model[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); @@ -1930,7 +1246,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, Reff = radi*radj/(radi+radj); bool touchflag; - if (normal == JKR){ + if (normal_model[itype][jtype] == JKR){ R2 = Reff*Reff; coh = normal_coeffs[itype][jtype][3]; a = cbrt(9.0*M_PI*coh*R2/(4*E)); @@ -2022,7 +1338,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, delta = radsum - r; dR = delta*Reff; - if (normal == JKR){ + if (normal_model[itype][jtype] == JKR){ dR2 = dR*dR; t0 = coh*coh*R2*R2*E; t1 = PI27SQ*t0; @@ -2042,22 +1358,21 @@ double PairGranular::single(int i, int j, int itype, int jtype, else{ knfac = E; Fne = knfac*delta; - if (normal != HOOKE) - a = sqrt(dR); + a = sqrt(dR); + if (normal_model[itype][jtype] != HOOKE) Fne *= a; - if (normal == DMT) + if (normal_model[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 == VELOCITY){ + if (damping_model[itype][jtype] == VELOCITY){ damp_normal = normal_coeffs[itype][jtype][1]; } - else if (damping == VISCOELASTIC){ - if (normal == HOOKE) a = sqrt(dR); + else if (damping_model[itype][jtype] == VISCOELASTIC){ damp_normal = normal_coeffs[itype][jtype][1]*a*meff; } - else if (damping == TSUJI){ + else if (damping_model[itype][jtype] == TSUJI){ damp_normal = normal_coeffs[itype][jtype][1]*sqrt(meff*knfac); } @@ -2099,17 +1414,19 @@ double PairGranular::single(int i, int j, int itype, int jtype, vrel = vtr1*vtr1 + vtr2*vtr2 + vtr3*vtr3; vrel = sqrt(vrel); - Fcrit = fabs(Fne); - if (normal == JKR){ + if (normal_model[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 != HOOKE){ + if (normal_model[itype][jtype] != HOOKE){ k_tangential *= a; } damp_tangential = tangential_coeffs[itype][jtype][1]*damp_normal; @@ -2150,7 +1467,7 @@ double PairGranular::single(int i, int j, int itype, int jtype, // Rolling resistance //**************************************** - if (roll != ROLL_NONE){ + if (roll_model[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]; @@ -2211,9 +1528,9 @@ double PairGranular::single(int i, int j, int itype, int jtype, //**************************************** // Twisting torque, including history effects //**************************************** - if (twist != TWIST_NONE){ + if (twist_model[itype][jtype] != TWIST_NONE){ magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of Marshall) - if (twist == TWIST_MARSHALL){ + if (twist_model[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; @@ -2293,22 +1610,18 @@ double PairGranular::memory_usage() mixing of Young's modulus (E) ------------------------------------------------------------------------- */ -double PairGranular::mix_stiffnessE(double Eii, double Ejj, double Gii, double Gjj) +double PairGranular::mix_stiffnessE(double Eii, double Ejj, double poisii, double poisjj) { - 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); + return 1/((1-poisii*poisii)/Eii+(1-poisjj*poisjj)/Ejj); } /* ---------------------------------------------------------------------- mixing of shear modulus (G) - ------------------------------------------------------------------------- */ +------------------------------------------------------------------------ */ -double PairGranular::mix_stiffnessG(double Eii, double Ejj, double Gii, double Gjj) +double PairGranular::mix_stiffnessG(double Gii, double Gjj, double poisii, double poisjj) { - 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); + return 1/((2.0 -poisii)/Gii+(2.0-poisjj)/Gjj); } /* ---------------------------------------------------------------------- @@ -2325,13 +1638,14 @@ double PairGranular::mix_geom(double valii, double valjj) Compute pull-off distance (beyond contact) for a given radius and atom type ------------------------------------------------------------------------- */ -double PairGranular::pulloff_distance(double radius, int itype) +double PairGranular::pulloff_distance(double radi, double radj, int itype, int jtype) { - double E, coh, a, delta_pulloff; + double E, coh, a, delta_pulloff, Reff; + Reff = radi*radj/(radi+radj); + if (Reff <= 0) return 0; 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); + E = normal_coeffs[itype][jtype][0]; + a = cbrt(9*M_PI*coh*Reff/(4*E)); + return a*a/Reff - 2*sqrt(M_PI*coh*a/E); } diff --git a/src/GRANULAR/pair_granular.h b/src/GRANULAR/pair_granular.h index f39f31e4cb..625ff17c72 100644 --- a/src/GRANULAR/pair_granular.h +++ b/src/GRANULAR/pair_granular.h @@ -21,40 +21,28 @@ PairStyle(granular,PairGranular) #define LMP_PAIR_GRANULAR_H #include "pair.h" +#include "pair_granular.h" namespace LAMMPS_NS { -class PairGranular : public Pair { +class PairGranular : public Pair{ public: PairGranular(class LAMMPS *); - virtual ~PairGranular(); - + ~PairGranular(); 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 settings(int, char **); + 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 &); + 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; @@ -72,22 +60,26 @@ public: double *mass_rigid; // rigid mass for owned+ghost atoms int nmax; // allocated size of mass_rigid - virtual void allocate(); + void allocate(); private: int size_history; - //Models - int normal, damping, tangential, roll, twist; + //Models choices + int **normal_model, **damping_model; + double **tangential_model, **roll_model, **twist_model; //History flags - int tangential_history, roll_history, twist_history; + int normal_history, tangential_history, roll_history, twist_history; //Indices of history entries - int tangential_history_index, roll_history_index, twist_history_index; + int normal_history_index; + int tangential_history_index; + int roll_history_index; + int twist_history_index; - //Flags for whether model choices have been set - int normal_set, tangential_set, damping_set, roll_set, twist_set; + //Per-type material coefficients + double **Emod, **poiss, **Gmod; //Per-type coefficients, set in pair coeff command double ***normal_coeffs; @@ -95,13 +87,14 @@ private: double ***roll_coeffs; double ***twist_coeffs; - //Optional user-specified global cutoff + //Optional user-specified global cutoff, per-type user-specified cutoffs + double **cutoff_type; double cutoff_global; - 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); + double mix_stiffnessE(double, double, double, double); + double mix_stiffnessG(double, double, double, double); + double mix_geom(double, double); + double pulloff_distance(double, double, int, int); }; } diff --git a/src/pair_granular.cpp b/src/pair_granular.cpp new file mode 100644 index 0000000000..82a470f83b --- /dev/null +++ b/src/pair_granular.cpp @@ -0,0 +1,2337 @@ +/* ---------------------------------------------------------------------- +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: Leo Silbert (SNL), Gary Grest (SNL), + Jeremy Lechman (SNL), Dan Bolintineanu (SNL), Ishan Srivastava (SNL) +----------------------------------------------------------------------- */ + +#include +#include +#include +#include +#include "pair_granular.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}; + +/* ---------------------------------------------------------------------- */ + +PairGranular::PairGranular(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; + +} + +/* ---------------------------------------------------------------------- */ +PairGranular::~PairGranular() +{ + 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); + + delete [] onerad_dynamic; + delete [] onerad_frozen; + delete [] maxrad_dynamic; + delete [] maxrad_frozen; + } + memory->destroy(mass_rigid); +} + +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_twist, int Tp_roll > +void PairGranular::compute_templated(int eflag, int vflag) +#else +void PairGranular::compute_untemplated + (int Tp_normal, int Tp_damping, int Tp_tangential, + int Tp_twist, int Tp_roll, 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; + + 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 (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); + } + 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 (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; + 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 (Tp_normal != HOOKE) + a = sqrt(dR); + Fne *= a; + 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 (Tp_damping == VELOCITY){ + damp_normal = 1; + } + else if (Tp_damping == VISCOELASTIC){ + if (Tp_normal == HOOKE) a = sqrt(dR); + damp_normal = a*meff; + } + else if (Tp_damping == TSUJI){ + damp_normal = sqrt(meff*knfac); + } + + Fdamp = -normal_coeffs[itype][jtype][1]*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 (Tp_normal == 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 (Tp_tangential > 0){ + 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 (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]; + + // 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 (Tp_roll > 1){ + 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 (Tp_twist != TWIST_NONE){ + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of 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; + } + else{ + k_twist = twist_coeffs[itype][jtype][0]; + damp_twist = twist_coeffs[itype][jtype][1]; + mu_twist = twist_coeffs[itype][jtype][2]; + } + if (Tp_twist > 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 (Tp_twist != TWIST_NONE){ + tortwist1 = magtortwist * nx; + tortwist2 = magtortwist * ny; + tortwist3 = magtortwist * nz; + + torque[i][0] += tortwist1; + torque[i][1] += tortwist2; + torque[i][2] += tortwist3; + } + + 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); + + 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 (Tp_twist != TWIST_NONE){ + torque[j][0] -= tortwist1; + torque[j][1] -= tortwist2; + torque[j][2] -= tortwist3; + } + if (Tp_roll != 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 PairGranular::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"); + + 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 PairGranular::settings(int narg, char **arg) +{ + if (narg == 1){ + cutoff_global = force->numeric(FLERR,arg[0]); + } + else{ + cutoff_global = -1; //Will be set based on particle sizes, model choice + } + tangential_history = 0; + roll_history = twist_history = 0; + normal_set = tangential_set = damping_set = roll_set = twist_set = 0; +} + +/* ---------------------------------------------------------------------- + set coeffs for one or more type pairs +------------------------------------------------------------------------- */ + +void PairGranular::coeff(int narg, char **arg) +{ + double normal_coeffs_local[4]; + double tangential_coeffs_local[4]; + double roll_coeffs_local[4]; + double twist_coeffs_local[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); + + 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"); + if (!normal_set) normal = HOOKE; + else if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); + 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"); + if (!normal_set) normal = HERTZ; + else if (normal_set && normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); + 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"); + if (!normal_set) normal = HERTZ_MATERIAL; + else if (normal != HERTZ) if (normal != HOOKE) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); + 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"); + if (!normal_set) normal = DMT; + else if (normal != DMT) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); + 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; + if (!normal_set) normal = JKR; + else if (normal != JKR) error->all(FLERR, "Illegal pair_coeff command, choice of normal contact model must be the same for all types"); + 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") == 0){ + if (iarg+1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters provided for damping model"); + if (strcmp(arg[iarg+1], "velocity") == 0){ + if (!damping_set) damping = VELOCITY; + else if (damping != VELOCITY) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "viscoelastic") == 0){ + if (!damping_set) damping = VISCOELASTIC; + else if (damping != VISCOELASTIC) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "tsuji") == 0){ + if (!damping_set) damping = TSUJI; + if (damping != TSUJI) error->all(FLERR, "Illegal pair_coeff command, choice of damping contact model must be the same for all types"); + } + damping_set = 1; + iarg += 2; + } + else if (strcmp(arg[iarg], "tangential") == 0){ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for tangential model"); + if (strcmp(arg[iarg+1], "nohistory") == 0){ + if (!tangential_set) tangential = TANGENTIAL_NOHISTORY; + else if (tangential != TANGENTIAL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "mindlin") == 0){ + if (!tangential_set) tangential = TANGENTIAL_MINDLIN; + else if (tangential != TANGENTIAL_MINDLIN) error->all(FLERR, "Illegal pair_coeff command, choice of tangential contact model must be the same for all types");; + 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. + iarg += 5; + } + else if (strcmp(arg[iarg], "rolling") == 0){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (strcmp(arg[iarg+1], "none") == 0){ + if (!roll_set) roll = ROLL_NONE; + else if (roll != ROLL_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); + iarg += 2; + } + else{ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for rolling model"); + if (strcmp(arg[iarg+1], "nohistory") == 0){ + if (!roll_set) roll = ROLL_NOHISTORY; + else if (roll != ROLL_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "sds") == 0){ + if (!roll_set) roll = ROLL_SDS; + else if (roll != ROLL_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of rolling friction model must be the same for all types"); + 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]); //kR + roll_coeffs_local[1] = force->numeric(FLERR,arg[iarg+3]); //gammaR + roll_coeffs_local[2] = force->numeric(FLERR,arg[iarg+4]); //rolling friction coeff. + iarg += 5; + } + } + else if (strcmp(arg[iarg], "twisting") == 0){ + if (iarg + 1 >= narg) error->all(FLERR, "Illegal pair_coeff command, not enough parameters"); + if (strcmp(arg[iarg+1], "none") == 0){ + if (!twist_set) twist = TWIST_NONE; + else if (twist != TWIST_NONE) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + iarg += 2; + } + else if (strcmp(arg[iarg+1], "marshall") == 0){ + if (!twist_set) twist = TWIST_MARSHALL; + else if (twist != TWIST_MARSHALL) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + twist_history = 1; + twist_set = 1; + iarg += 2; + } + else{ + if (iarg + 4 >= narg) error->all(FLERR,"Illegal pair_coeff command, not enough parameters provided for twist model"); + else if (strcmp(arg[iarg+1], "nohistory") == 0){ + if (!twist_set) twist = TWIST_NOHISTORY; + if (twist != TWIST_NOHISTORY) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + } + else if (strcmp(arg[iarg+1], "sds") == 0){ + if (!twist_set) twist = TWIST_SDS; + else if (twist != TWIST_SDS) error->all(FLERR, "Illegal pair_coeff command, choice of twisting friction model must be the same for all types"); + 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. + iarg += 5; + } + } + else error->all(FLERR, "Illegal pair coeff command"); + } + + //It is an error not to specify normal or tangential model + if (!normal_set) error->all(FLERR, "Illegal pair coeff command, must specify normal contact model"); + if (!tangential_set) error->all(FLERR, "Illegal pair coeff command, must specify tangential contact model"); + + //If unspecified, set damping to VISCOELASTIC, twist/roll to NONE (cannot be changed by subsequent pair_coeff commands) + if (!damping_set) damping = VISCOELASTIC; + if (!roll_set) roll = ROLL_NONE; + if (!twist_set) twist = TWIST_NONE; + damping_set = roll_set = twist_set = 1; + + int count = 0; + double damp; + 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++) { + normal_coeffs[i][j][0] = normal_coeffs[j][i][0] = normal_coeffs_local[0]; + normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = damp; + 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[j][i][3] = normal_coeffs_local[3]; + + for (int k = 0; k < 3; k++) + tangential_coeffs[i][j][k] = tangential_coeffs[j][i][k] = tangential_coeffs_local[k]; + + if (roll != ROLL_NONE) + for (int k = 0; k < 3; k++) + roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = roll_coeffs_local[k]; + + if (twist != TWIST_NONE && twist != TWIST_MARSHALL) + for (int k = 0; k < 3; k++) + twist_coeffs[i][j][k] = twist_coeffs[j][i][k] = twist_coeffs_local[k]; + + setflag[i][j] = 1; + count++; + } + } + + if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients"); +} + +/* ---------------------------------------------------------------------- + init specific to this pair style +------------------------------------------------------------------------- */ + +void PairGranular::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 + if (normal == 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 == 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 == 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 == 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 PairGranular::init_one(int i, int j) +{ + double cutoff; + if (setflag[i][j] == 0) { + + if (normal != HOOKE && normal != HERTZ){ + normal_coeffs[i][j][0] = normal_coeffs[j][i][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] = normal_coeffs[j][i][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] = normal_coeffs[j][i][0] = mix_geom(normal_coeffs[i][i][0], normal_coeffs[j][j][0]); + } + + normal_coeffs[i][j][1] = normal_coeffs[j][i][1] = mix_geom(normal_coeffs[i][i][1], normal_coeffs[j][j][1]); + if ((normal == JKR) || (normal == DMT)) + normal_coeffs[i][j][3] = normal_coeffs[j][i][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] = tangential_coeffs[j][i][k] = mix_geom(tangential_coeffs[i][i][k], tangential_coeffs[j][j][k]); + + + if (roll != ROLL_NONE){ + for (int k = 0; k < 3; k++) + roll_coeffs[i][j][k] = roll_coeffs[j][i][k] = mix_geom(roll_coeffs[i][i][k], roll_coeffs[j][j][k]); + } + + if (twist != TWIST_NONE && twist != TWIST_MARSHALL){ + for (int k = 0; k < 3; k++) + twist_coeffs[i][j][k] = twist_coeffs[j][i][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 (cutoff_global < 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; + } + } + else{ + cutoff = cutoff_global; + } + return cutoff; +} + + +/* ---------------------------------------------------------------------- + proc 0 writes to restart file + ------------------------------------------------------------------------- */ + +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_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 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_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_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 PairGranular::reset_dt() +{ + dt = update->dt; +} + +/* ---------------------------------------------------------------------- */ + +double PairGranular::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 == 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 == 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 != HOOKE) + a = sqrt(dR); + Fne *= a; + 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 == VELOCITY){ + damp_normal = normal_coeffs[itype][jtype][1]; + } + else if (damping == VISCOELASTIC){ + if (normal == HOOKE) a = sqrt(dR); + damp_normal = normal_coeffs[itype][jtype][1]*a*meff; + } + else if (damping == 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); + + Fcrit = fabs(Fne); + if (normal == JKR){ + F_pulloff = 3*M_PI*coh*Reff; + Fcrit = fabs(Fne + 2*F_pulloff); + } + + //------------------------------ + //Tangential forces + //------------------------------ + k_tangential = tangential_coeffs[itype][jtype][0]; + if (normal != 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 != 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 != TWIST_NONE){ + magtwist = relrot1*nx + relrot2*ny + relrot3*nz; //Omega_T (eq 29 of 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; + } + 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 PairGranular::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 PairGranular::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 PairGranular::memory_usage() +{ + double bytes = nmax * sizeof(double); + return bytes; +} + +/* ---------------------------------------------------------------------- + mixing of Young's modulus (E) +------------------------------------------------------------------------- */ + +double PairGranular::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 PairGranular::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 PairGranular::mix_geom(double valii, double valjj) +{ + return sqrt(valii*valjj); +} + + +/* ---------------------------------------------------------------------- + Compute pull-off distance (beyond contact) for a given radius and atom type +------------------------------------------------------------------------- */ + +double PairGranular::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/pair_granular.h b/src/pair_granular.h new file mode 100644 index 0000000000..f39f31e4cb --- /dev/null +++ b/src/pair_granular.h @@ -0,0 +1,120 @@ +/* ---------------------------------------------------------- + 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,PairGranular) + +#else + +#ifndef LMP_PAIR_GRANULAR_H +#define LMP_PAIR_GRANULAR_H + +#include "pair.h" + +namespace LAMMPS_NS { + +class PairGranular : public Pair { +public: + PairGranular(class LAMMPS *); + virtual ~PairGranular(); + + 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(); + 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; + + //Models + int normal, damping, tangential, roll, twist; + + //History flags + int tangential_history, roll_history, twist_history; + + //Indices of history entries + int tangential_history_index, roll_history_index, twist_history_index; + + //Flags for whether model choices have been set + int normal_set, tangential_set, damping_set, roll_set, twist_set; + + //Per-type coefficients, set in pair coeff command + double ***normal_coeffs; + double ***tangential_coeffs; + double ***roll_coeffs; + double ***twist_coeffs; + + //Optional user-specified global cutoff + double cutoff_global; + + 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. + + */