1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
/* -*- c++ -*- ----------------------------------------------------------
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(body/rounded/polyhedron,PairBodyRoundedPolyhedron)
#else
#ifndef LMP_PAIR_BODY_ROUNDED_POLYHEDRON_H
#define LMP_PAIR_BODY_ROUNDED_POLYHEDRON_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBodyRoundedPolyhedron : public Pair {
public:
PairBodyRoundedPolyhedron(class LAMMPS *);
~PairBodyRoundedPolyhedron();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
virtual void kernel_force(double R, int itype, int jtype,
double& energy, double& fpair);
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
int type; // 0 = VERTEX-FACE; 1 = EDGE-EDGE
double fx,fy,fz; // unscaled cohesive forces at contact
double xi[3]; // coordinates of the contact point on ibody
double xj[3]; // coordinates of the contact point on jbody
double separation; // contact surface separation
int unique;
};
protected:
double **k_n; // normal repulsion strength
double **k_na; // normal attraction strength
double c_n; // normal damping coefficient
double c_t; // tangential damping coefficient
double mu; // normal friction coefficient during gross sliding
double A_ua; // characteristic contact area
double cut_inner; // cutoff for interaction between vertex-edge surfaces
class AtomVecBody *avec;
class BodyRoundedPolyhedron *bptr;
double **discrete; // list of all sub-particles for all bodies
int ndiscrete; // number of discretes in list
int dmax; // allocated size of discrete list
int *dnum; // number of discretes per line, 0 if uninit
int *dfirst; // index of first discrete per each line
int nmax; // allocated size of dnum,dfirst vectors
double **edge; // list of all edge for all bodies
int nedge; // number of edge in list
int edmax; // allocated size of edge list
int *ednum; // number of edges per line, 0 if uninit
int *edfirst; // index of first edge per each line
int ednummax; // allocated size of ednum,edfirst vectors
double **face; // list of all edge for all bodies
int nface; // number of faces in list
int facmax; // allocated size of face list
int *facnum; // number of faces per line, 0 if uninit
int *facfirst; // index of first face per each line
int facnummax; // allocated size of facnum,facfirst vectors
double *enclosing_radius; // enclosing radii for all bodies
double *rounded_radius; // rounded radii for all bodies
double *maxerad; // per-type maximum enclosing radius
void allocate();
void body2space(int);
// sphere-sphere interaction
void sphere_against_sphere(int ibody, int jbody, int itype, int jtype,
double delx, double dely, double delz, double rsq,
double** v, double** f, int evflag);
// sphere-edge interaction
void sphere_against_edge(int ibody, int jbody, int itype, int jtype,
double** x, double** v, double** f, double** torque,
double** angmom, int evflag);
// sphere-face interaction
void sphere_against_face(int ibody, int jbody, int itype, int jtype,
double** x, double** v, double** f, double** torque,
double** angmom, int evflag);
// edge-edge interactions
int edge_against_edge(int ibody, int jbody, int itype, int jtype,
double** x,Contact* contact_list, int &num_contacts,
double &evdwl, double* facc);
// edge-face interactions
int edge_against_face(int ibody, int jbody, int itype, int jtype,
double** x, Contact* contact_list, int &num_contacts,
double &evdwl, double* facc);
// a face vs. a single edge
int interaction_face_to_edge(int ibody, int face_index, double* xmi,
double rounded_radius_i, int jbody, int edge_index,
double* xmj, double rounded_radius_j,
int itype, int jtype, double cut_inner,
Contact* contact_list, int &num_contacts,
double& energy, double* facc);
// an edge vs. an edge from another body
int interaction_edge_to_edge(int ibody, int edge_index_i, double* xmi,
double rounded_radius_i, int jbody, int edge_index_j,
double* xmj, double rounded_radius_j,
int itype, int jtype, double cut_inner,
Contact* contact_list, int &num_contacts,
double& energy, double* facc);
// compute contact forces if contact points are detected
void contact_forces(int ibody, int jbody, double *xi, double *xj,
double delx, double dely, double delz, double fx, double fy, double fz,
double** x, double** v, double** angmom, double** f, double** torque,
double* facc);
// compute force and torque between two bodies given a pair of interacting points
void pair_force_and_torque(int ibody, int jbody, double* pi, double* pj,
double r, double contact_dist, int itype, int jtype,
double** x, double** v, double** f, double** torque,
double** angmom, int jflag, double& energy, double* facc);
// rescale the cohesive forces if a contact area is detected
void rescale_cohesive_forces(double** x, double** f, double** torque,
Contact* contact_list, int &num_contacts,
int itype, int jtype, double* facc);
// compute the separation between two contacts
double contact_separation(const Contact& c1, const Contact& c2);
// detect the unique contact points (as there may be double counts)
void find_unique_contacts(Contact* contact_list, int& num_contacts);
// accumulate torque to a body given a force at a given point
void sum_torque(double* xm, double *x, double fx, double fy, double fz, double* torque);
// find the intersection point (if any) between an edge and a face
int edge_face_intersect(double* x1, double* x2, double* x3, double* a, double* b,
double* hi1, double* hi2, double& d1, double& d2,
int& inside_a, int& inside_b);
// helper functions
int opposite_sides(double* n, double* x0, double* a, double* b);
void project_pt_plane(const double* q, const double* p,
const double* n, double* q_proj, double &d);
void project_pt_plane(const double* q, const double* x1, const double* x2,
const double* x3, double* q_proj, double &d, int& inside);
void project_pt_line(const double* q, const double* xi1, const double* xi2,
double* h, double& d, double& t);
void inside_polygon(int ibody, int face_index, double* xmi,
const double* q1, const double* q2, int& inside1, int& inside2);
void distance_bt_edges(const double* x1, const double* x2,
const double* x3, const double* x4,
double* h1, double* h2, double& t1, double& t2, double& r);
void total_velocity(double* p, double *xcm, double* vcm, double *angmom,
double *inertia, double *quat, double* vi);
void sanity_check();
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair body/rounded/polyhedron requires atom style body rounded/polyhedron
Self-explanatory.
E: Pair body requires body style rounded/polyhedron
This pair style is specific to the rounded/polyhedron body style.
*/