1 // clang-format off
2 /* ----------------------------------------------------------------------
3    LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
4    https://www.lammps.org/, Sandia National Laboratories
5    Steve Plimpton, sjplimp@sandia.gov
6 
7    Copyright (2003) Sandia Corporation.  Under the terms of Contract
8    DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9    certain rights in this software.  This software is distributed under
10    the GNU General Public License.
11 
12    See the README file in the top-level LAMMPS directory.
13 ------------------------------------------------------------------------- */
14 
15 /* ----------------------------------------------------------------------
16    Contributing author: Trung Dac Nguyen (ndactrung@gmail.com)
17    Ref: Fraige, Langston, Matchett and Dodds, Particuology 2008, 6:455-466
18    Note: The current implementation has not taken into account
19          the contact history for friction forces.
20 ------------------------------------------------------------------------- */
21 
22 #include "pair_body_rounded_polygon.h"
23 
24 #include "atom.h"
25 #include "atom_vec_body.h"
26 #include "body_rounded_polygon.h"
27 #include "comm.h"
28 #include "error.h"
29 #include "fix.h"
30 #include "force.h"
31 #include "math_extra.h"
32 #include "memory.h"
33 #include "modify.h"
34 #include "neigh_list.h"
35 #include "neighbor.h"
36 
37 #include <cmath>
38 #include <cstring>
39 
40 using namespace LAMMPS_NS;
41 
42 #define DELTA 10000
43 #define EPSILON 1e-3
44 #define MAX_CONTACTS 4  // maximum number of contacts for 2D models
45 #define EFF_CONTACTS 2  // effective contacts for 2D models
46 
47 //#define _CONVEX_POLYGON
48 //#define _POLYGON_DEBUG
49 
50 enum {INVALID=0,NONE=1,VERTEXI=2,VERTEXJ=3,EDGE=4};
51 
52 /* ---------------------------------------------------------------------- */
53 
PairBodyRoundedPolygon(LAMMPS * lmp)54 PairBodyRoundedPolygon::PairBodyRoundedPolygon(LAMMPS *lmp) : Pair(lmp)
55 {
56   dmax = nmax = 0;
57   discrete = nullptr;
58   dnum = dfirst = nullptr;
59 
60   edmax = ednummax = 0;
61   edge = nullptr;
62   ednum = edfirst = nullptr;
63 
64   enclosing_radius = nullptr;
65   rounded_radius = nullptr;
66   maxerad = nullptr;
67 
68   single_enable = 0;
69   restartinfo = 0;
70 
71   c_n = 0.1;
72   c_t = 0.2;
73   mu = 0.0;
74   delta_ua = 1.0;
75 }
76 
77 /* ---------------------------------------------------------------------- */
78 
~PairBodyRoundedPolygon()79 PairBodyRoundedPolygon::~PairBodyRoundedPolygon()
80 {
81   memory->destroy(discrete);
82   memory->destroy(dnum);
83   memory->destroy(dfirst);
84 
85   memory->destroy(edge);
86   memory->destroy(ednum);
87   memory->destroy(edfirst);
88 
89   memory->destroy(enclosing_radius);
90   memory->destroy(rounded_radius);
91 
92   if (allocated) {
93     memory->destroy(setflag);
94     memory->destroy(cutsq);
95 
96     memory->destroy(k_n);
97     memory->destroy(k_na);
98     memory->destroy(maxerad);
99   }
100 }
101 
102 /* ---------------------------------------------------------------------- */
103 
compute(int eflag,int vflag)104 void PairBodyRoundedPolygon::compute(int eflag, int vflag)
105 {
106   int i,j,ii,jj,inum,jnum,itype,jtype;
107   int ni,nj,npi,npj,ifirst,jfirst;
108   int nei,nej,iefirst,jefirst;
109   double xtmp,ytmp,ztmp,delx,dely,delz,evdwl;
110   double rsq,r,radi,radj,k_nij,k_naij;
111   double facc[3];
112   int *ilist,*jlist,*numneigh,**firstneigh;
113 
114   evdwl = 0.0;
115   ev_init(eflag,vflag);
116 
117   double **x = atom->x;
118   double **v = atom->v;
119   double **f = atom->f;
120   double **torque = atom->torque;
121   double **angmom = atom->angmom;
122   double *radius = atom->radius;
123   tagint* tag = atom->tag;
124   int *body = atom->body;
125   int *type = atom->type;
126   int nlocal = atom->nlocal;
127   int nall = nlocal + atom->nghost;
128   int newton_pair = force->newton_pair;
129 
130   inum = list->inum;
131   ilist = list->ilist;
132   numneigh = list->numneigh;
133   firstneigh = list->firstneigh;
134 
135   // grow the per-atom lists if necessary and initialize
136 
137   if (atom->nmax > nmax) {
138     memory->destroy(dnum);
139     memory->destroy(dfirst);
140     memory->destroy(ednum);
141     memory->destroy(edfirst);
142     memory->destroy(enclosing_radius);
143     memory->destroy(rounded_radius);
144     nmax = atom->nmax;
145     memory->create(dnum,nmax,"pair:dnum");
146     memory->create(dfirst,nmax,"pair:dfirst");
147     memory->create(ednum,nmax,"pair:ednum");
148     memory->create(edfirst,nmax,"pair:edfirst");
149     memory->create(enclosing_radius,nmax,"pair:enclosing_radius");
150     memory->create(rounded_radius,nmax,"pair:rounded_radius");
151   }
152 
153   ndiscrete = nedge = 0;
154   for (i = 0; i < nall; i++)
155     dnum[i] = ednum[i] = 0;
156 
157   // loop over neighbors of my atoms
158 
159   for (ii = 0; ii < inum; ii++) {
160     i = ilist[ii];
161     xtmp = x[i][0];
162     ytmp = x[i][1];
163     ztmp = x[i][2];
164     itype = type[i];
165     radi = radius[i];
166     jlist = firstneigh[i];
167     jnum = numneigh[i];
168 
169     if (body[i] >= 0) {
170       if (dnum[i] == 0) body2space(i);
171       npi = dnum[i];
172       ifirst = dfirst[i];
173       nei = ednum[i];
174       iefirst = edfirst[i];
175     }
176 
177     for (jj = 0; jj < jnum; jj++) {
178       j = jlist[jj];
179       j &= NEIGHMASK;
180 
181       delx = xtmp - x[j][0];
182       dely = ytmp - x[j][1];
183       delz = ztmp - x[j][2];
184       rsq = delx*delx + dely*dely + delz*delz;
185       jtype = type[j];
186       radj = radius[j];
187 
188       // body/body interactions
189 
190       evdwl = 0.0;
191       facc[0] = facc[1] = facc[2] = 0;
192 
193       if (body[i] < 0 || body[j] < 0) continue;
194 
195       if (dnum[j] == 0) body2space(j);
196       npj = dnum[j];
197       jfirst = dfirst[j];
198       nej = ednum[j];
199       jefirst = edfirst[j];
200 
201       k_nij = k_n[itype][jtype];
202       k_naij = k_na[itype][jtype];
203 
204       // no interaction
205 
206       r = sqrt(rsq);
207       if (r > radi + radj + cut_inner) continue;
208 
209       if (npi == 1 && npj == 1) {
210         sphere_against_sphere(i, j, delx, dely, delz, rsq, k_nij, k_naij, x, v, f, evflag);
211         continue;
212       }
213 
214       // reset vertex and edge forces
215 
216       for (ni = 0; ni < npi; ni++) {
217         discrete[ifirst+ni][3] = 0;
218         discrete[ifirst+ni][4] = 0;
219         discrete[ifirst+ni][5] = 0;
220       }
221 
222       for (nj = 0; nj < npj; nj++) {
223         discrete[jfirst+nj][3] = 0;
224         discrete[jfirst+nj][4] = 0;
225         discrete[jfirst+nj][5] = 0;
226       }
227 
228       for (ni = 0; ni < nei; ni++) {
229         edge[iefirst+ni][2] = 0;
230         edge[iefirst+ni][3] = 0;
231         edge[iefirst+ni][4] = 0;
232       }
233 
234       for (nj = 0; nj < nej; nj++) {
235         edge[jefirst+nj][2] = 0;
236         edge[jefirst+nj][3] = 0;
237         edge[jefirst+nj][4] = 0;
238       }
239 
240       int num_contacts, done;
241       double delta_a, j_a;
242       Contact contact_list[MAX_CONTACTS];
243 
244       num_contacts = 0;
245 
246       // check interaction between i's vertices and j' edges
247 
248       vertex_against_edge(i, j, k_nij, k_naij, x, f, torque, tag,
249                           contact_list, num_contacts, evdwl, facc);
250 
251       // check interaction between j's vertices and i' edges
252 
253       vertex_against_edge(j, i, k_nij, k_naij, x, f, torque, tag,
254                           contact_list, num_contacts, evdwl, facc);
255 
256       if (num_contacts >= 2) {
257 
258         // find the first two distinct contacts
259 
260         done = 0;
261         for (int m = 0; m < num_contacts-1; m++) {
262           for (int n = m+1; n < num_contacts; n++) {
263             delta_a = contact_separation(contact_list[m], contact_list[n]);
264             if (delta_a > 0) {
265               j_a = delta_a / (EFF_CONTACTS * delta_ua);
266               if (j_a < 1.0) j_a = 1.0;
267 
268               // scale the force at both contacts
269 
270               contact_forces(contact_list[m], j_a, x, v, angmom, f, torque,
271                              evdwl, facc);
272               contact_forces(contact_list[n], j_a, x, v, angmom, f, torque,
273                              evdwl, facc);
274               done = 1;
275 
276               #ifdef _POLYGON_DEBUG
277               printf("  Two separate contacts %d and %d: delta_a = %f; j_a = %f\n",
278                 m, n, delta_a, j_a);
279               printf("    %d: vertex %d of body %d and edge %d of body %d; "
280                      "xv = %f %f %f; xe = %f %f %f\n",
281                      m, contact_list[m].vertex, contact_list[m].ibody,
282                      contact_list[m].edge, contact_list[m].jbody,
283                      contact_list[m].xv[0], contact_list[m].xv[1],
284                      contact_list[m].xv[2], contact_list[m].xe[0],
285                      contact_list[m].xe[1], contact_list[m].xe[2]);
286               printf("    %d: vertex %d of body %d and edge %d of body %d; "
287                      "xv = %f %f %f; xe = %f %f %f\n",
288                      n, contact_list[n].vertex, contact_list[n].ibody,
289                      contact_list[n].edge, contact_list[n].jbody,
290                      contact_list[n].xv[0], contact_list[n].xv[1],
291                      contact_list[n].xv[2], contact_list[n].xe[0],
292                      contact_list[n].xe[1], contact_list[n].xe[2]);
293               #endif
294 
295               break;
296             }
297           }
298           if (done == 1) break;
299         }
300 
301 
302       } else if (num_contacts == 1) {
303 
304         // if there's only one contact, it should be handled here
305         // since forces/torques have not been accumulated from vertex2edge()
306 
307         contact_forces(contact_list[0], 1.0, x, v, angmom, f, torque, evdwl, facc);
308 
309         #ifdef _POLYGON_DEBUG
310         printf("One contact between vertex %d of body %d and edge %d of body %d:\n",
311                 contact_list[0].vertex, tag[contact_list[0].ibody],
312                 contact_list[0].edge, tag[contact_list[0].jbody]);
313         printf("xv = %f %f %f; xe = %f %f %f\n",
314                contact_list[0].xv[0], contact_list[0].xv[1], contact_list[0].xv[2],
315                contact_list[0].xe[0], contact_list[0].xe[1], contact_list[0].xe[2]);
316         #endif
317       }
318 
319       #ifdef _POLYGON_DEBUG
320       int num_overlapping_contacts = 0;
321       for (int m = 0; m < num_contacts-1; m++) {
322         for (int n = m+1; n < num_contacts; n++) {
323           double l = contact_separation(contact_list[m], contact_list[n]);
324           if (l < EPSILON) num_overlapping_contacts++;
325         }
326       }
327       printf("There are %d contacts detected, %d of which overlap.\n",
328              num_contacts, num_overlapping_contacts);
329       #endif
330 
331       if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,evdwl,0.0,
332                                facc[0],facc[1],facc[2],delx,dely,delz);
333 
334     } // end for jj
335   }
336 
337   if (vflag_fdotr) virial_fdotr_compute();
338 }
339 
340 /* ----------------------------------------------------------------------
341    allocate all arrays
342 ------------------------------------------------------------------------- */
343 
allocate()344 void PairBodyRoundedPolygon::allocate()
345 {
346   allocated = 1;
347   int n = atom->ntypes;
348 
349   memory->create(setflag,n+1,n+1,"pair:setflag");
350   for (int i = 1; i <= n; i++)
351     for (int j = i; j <= n; j++)
352       setflag[i][j] = 0;
353 
354   memory->create(cutsq,n+1,n+1,"pair:cutsq");
355 
356   memory->create(k_n,n+1,n+1,"pair:k_n");
357   memory->create(k_na,n+1,n+1,"pair:k_na");
358   memory->create(maxerad,n+1,"pair:maxerad");
359 }
360 
361 /* ----------------------------------------------------------------------
362    global settings
363 ------------------------------------------------------------------------- */
364 
settings(int narg,char ** arg)365 void PairBodyRoundedPolygon::settings(int narg, char **arg)
366 {
367   if (narg < 5) error->all(FLERR,"Illegal pair_style command");
368 
369   c_n = utils::numeric(FLERR,arg[0],false,lmp);
370   c_t = utils::numeric(FLERR,arg[1],false,lmp);
371   mu = utils::numeric(FLERR,arg[2],false,lmp);
372   delta_ua = utils::numeric(FLERR,arg[3],false,lmp);
373   cut_inner = utils::numeric(FLERR,arg[4],false,lmp);
374 
375   if (delta_ua < 0) delta_ua = 1;
376 }
377 
378 /* ----------------------------------------------------------------------
379    set coeffs for one or more type pairs
380 ------------------------------------------------------------------------- */
381 
coeff(int narg,char ** arg)382 void PairBodyRoundedPolygon::coeff(int narg, char **arg)
383 {
384   if (narg < 4 || narg > 5)
385     error->all(FLERR,"Incorrect args for pair coefficients");
386   if (!allocated) allocate();
387 
388   int ilo,ihi,jlo,jhi;
389   utils::bounds(FLERR,arg[0],1,atom->ntypes,ilo,ihi,error);
390   utils::bounds(FLERR,arg[1],1,atom->ntypes,jlo,jhi,error);
391 
392   double k_n_one = utils::numeric(FLERR,arg[2],false,lmp);
393   double k_na_one = utils::numeric(FLERR,arg[3],false,lmp);
394 
395   int count = 0;
396   for (int i = ilo; i <= ihi; i++) {
397     for (int j = MAX(jlo,i); j <= jhi; j++) {
398       k_n[i][j] = k_n_one;
399       k_na[i][j] = k_na_one;
400       setflag[i][j] = 1;
401       count++;
402     }
403   }
404 
405   if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
406 }
407 
408 /* ----------------------------------------------------------------------
409    init specific to this pair style
410 ------------------------------------------------------------------------- */
411 
init_style()412 void PairBodyRoundedPolygon::init_style()
413 {
414   avec = (AtomVecBody *) atom->style_match("body");
415   if (!avec)
416     error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
417   if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
418     error->all(FLERR,"Pair body/rounded/polygon requires "
419                "body style rounded/polygon");
420   bptr = (BodyRoundedPolygon *) avec->bptr;
421 
422   if (force->newton_pair == 0)
423     error->all(FLERR,"Pair style body/rounded/polygon requires "
424                "newton pair on");
425 
426   if (comm->ghost_velocity == 0)
427     error->all(FLERR,"Pair body/rounded/polygon requires "
428                "ghost atoms store velocity");
429 
430   neighbor->request(this);
431 
432   // find the maximum enclosing radius for each atom type
433 
434   int i, itype;
435   double eradi;
436   int* body = atom->body;
437   int* type = atom->type;
438   int ntypes = atom->ntypes;
439   int nlocal = atom->nlocal;
440 
441   if (atom->nmax > nmax) {
442     memory->destroy(dnum);
443     memory->destroy(dfirst);
444     memory->destroy(ednum);
445     memory->destroy(edfirst);
446     memory->destroy(enclosing_radius);
447     memory->destroy(rounded_radius);
448     nmax = atom->nmax;
449     memory->create(dnum,nmax,"pair:dnum");
450     memory->create(dfirst,nmax,"pair:dfirst");
451     memory->create(ednum,nmax,"pair:ednum");
452     memory->create(edfirst,nmax,"pair:edfirst");
453     memory->create(enclosing_radius,nmax,"pair:enclosing_radius");
454     memory->create(rounded_radius,nmax,"pair:rounded_radius");
455   }
456 
457   ndiscrete = nedge = 0;
458   for (i = 0; i < nlocal; i++)
459     dnum[i] = ednum[i] = 0;
460 
461   double *merad = nullptr;
462   memory->create(merad,ntypes+1,"pair:merad");
463   for (i = 1; i <= ntypes; i++)
464     maxerad[i] = merad[i] = 0;
465 
466   int ipour;
467   for (ipour = 0; ipour < modify->nfix; ipour++)
468     if (strcmp(modify->fix[ipour]->style,"pour") == 0) break;
469   if (ipour == modify->nfix) ipour = -1;
470 
471   int idep;
472   for (idep = 0; idep < modify->nfix; idep++)
473     if (strcmp(modify->fix[idep]->style,"deposit") == 0) break;
474   if (idep == modify->nfix) idep = -1;
475 
476   for (i = 1; i <= ntypes; i++) {
477     merad[i] = 0.0;
478     if (ipour >= 0) {
479       itype = i;
480       merad[i] =
481         *((double *) modify->fix[ipour]->extract("radius",itype));
482     }
483     if (idep >= 0) {
484       itype = i;
485       merad[i] =
486         *((double *) modify->fix[idep]->extract("radius",itype));
487     }
488   }
489 
490   for (i = 0; i < nlocal; i++) {
491     itype = type[i];
492     if (body[i] >= 0) {
493       if (dnum[i] == 0) body2space(i);
494       eradi = enclosing_radius[i];
495       if (eradi > merad[itype]) merad[itype] = eradi;
496     } else
497       merad[itype] = 0;
498   }
499 
500   MPI_Allreduce(&merad[1],&maxerad[1],ntypes,MPI_DOUBLE,MPI_MAX,world);
501 
502   memory->destroy(merad);
503 }
504 
505 /* ----------------------------------------------------------------------
506    init for one type pair i,j and corresponding j,i
507 ------------------------------------------------------------------------- */
508 
init_one(int i,int j)509 double PairBodyRoundedPolygon::init_one(int i, int j)
510 {
511   k_n[j][i] = k_n[i][j];
512   k_na[j][i] = k_na[i][j];
513 
514   return (maxerad[i]+maxerad[j]);
515 }
516 
517 /* ----------------------------------------------------------------------
518    convert N sub-particles in body I to space frame using current quaternion
519    store sub-particle space-frame displacements from COM in discrete list
520 ------------------------------------------------------------------------- */
521 
body2space(int i)522 void PairBodyRoundedPolygon::body2space(int i)
523 {
524   int ibonus = atom->body[i];
525   AtomVecBody::Bonus *bonus = &avec->bonus[ibonus];
526   int nsub = bptr->nsub(bonus);
527   double *coords = bptr->coords(bonus);
528   int body_num_edges = bptr->nedges(bonus);
529   double* edge_ends = bptr->edges(bonus);
530   double eradius = bptr->enclosing_radius(bonus);
531   double rradius = bptr->rounded_radius(bonus);
532 
533   // get the number of sub-particles (vertices)
534   // and the index of the first vertex of my body in the list
535 
536   dnum[i] = nsub;
537   dfirst[i] = ndiscrete;
538 
539   // grow the vertex list if necessary
540   // the first 3 columns are for coords, the last 3 for forces
541 
542   if (ndiscrete + nsub > dmax) {
543     dmax += DELTA;
544     memory->grow(discrete,dmax,6,"pair:discrete");
545   }
546 
547   double p[3][3];
548   MathExtra::quat_to_mat(bonus->quat,p);
549 
550   for (int m = 0; m < nsub; m++) {
551     MathExtra::matvec(p,&coords[3*m],discrete[ndiscrete]);
552     discrete[ndiscrete][3] = 0;
553     discrete[ndiscrete][4] = 0;
554     discrete[ndiscrete][5] = 0;
555     ndiscrete++;
556   }
557 
558   // get the number of edges (vertices)
559   // and the index of the first edge of my body in the list
560 
561   ednum[i] = body_num_edges;
562   edfirst[i] = nedge;
563 
564   // grow the edge list if necessary
565   // the first 2 columns are for vertex indices within body, the last 3 for forces
566 
567   if (nedge + body_num_edges > edmax) {
568     edmax += DELTA;
569     memory->grow(edge,edmax,5,"pair:edge");
570   }
571 
572   if ((body_num_edges > 0) && (edge_ends == nullptr))
573     error->one(FLERR,"Inconsistent edge data for body of atom {}",
574                                  atom->tag[i]);
575 
576   for (int m = 0; m < body_num_edges; m++) {
577     edge[nedge][0] = static_cast<int>(edge_ends[2*m+0]);
578     edge[nedge][1] = static_cast<int>(edge_ends[2*m+1]);
579     edge[nedge][2] = 0;
580     edge[nedge][3] = 0;
581     edge[nedge][4] = 0;
582     nedge++;
583   }
584 
585   enclosing_radius[i] = eradius;
586   rounded_radius[i] = rradius;
587 }
588 
589 /* ----------------------------------------------------------------------
590    Interaction between two spheres with different radii
591    according to the 2D model from Fraige et al.
592 ---------------------------------------------------------------------- */
593 
sphere_against_sphere(int i,int j,double delx,double dely,double delz,double rsq,double k_n,double k_na,double **,double ** v,double ** f,int evflag)594 void PairBodyRoundedPolygon::sphere_against_sphere(int i, int j,
595                        double delx, double dely, double delz, double rsq,
596                        double k_n, double k_na, double** /*x*/, double** v,
597                        double** f, int evflag)
598 {
599   double rradi,rradj;
600   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
601   double rij,rsqinv,R,fx,fy,fz,fn[3],ft[3],fpair,shift,energy;
602   int nlocal = atom->nlocal;
603   int newton_pair = force->newton_pair;
604 
605   rradi = rounded_radius[i];
606   rradj = rounded_radius[j];
607 
608   rsqinv = 1.0/rsq;
609   rij = sqrt(rsq);
610   R = rij - (rradi + rradj);
611   shift = k_na * cut_inner;
612 
613   energy = 0;
614 
615   if (R <= 0) {           // deformation occurs
616     fpair = -k_n * R - shift;
617     energy = (0.5 * k_n * R + shift) * R;
618   } else if (R <= cut_inner) {   // not deforming but cohesive ranges overlap
619     fpair = k_na * R - shift;
620     energy = (-0.5 * k_na * R + shift) * R;
621   } else fpair = 0.0;
622 
623   fx = delx*fpair/rij;
624   fy = dely*fpair/rij;
625   fz = delz*fpair/rij;
626 
627   if (R <= EPSILON) { // in contact
628 
629     // relative translational velocity
630 
631     vr1 = v[i][0] - v[j][0];
632     vr2 = v[i][1] - v[j][1];
633     vr3 = v[i][2] - v[j][2];
634 
635     // normal component
636 
637     vnnr = vr1*delx + vr2*dely + vr3*delz;
638     vn1 = delx*vnnr * rsqinv;
639     vn2 = dely*vnnr * rsqinv;
640     vn3 = delz*vnnr * rsqinv;
641 
642     // tangential component
643 
644     vt1 = vr1 - vn1;
645     vt2 = vr2 - vn2;
646     vt3 = vr3 - vn3;
647 
648     // normal friction term at contact
649 
650     fn[0] = -c_n * vn1;
651     fn[1] = -c_n * vn2;
652     fn[2] = -c_n * vn3;
653 
654     // tangential friction term at contact,
655     // excluding the tangential deformation term for now
656 
657     ft[0] = -c_t * vt1;
658     ft[1] = -c_t * vt2;
659     ft[2] = -c_t * vt3;
660 
661     fx += fn[0] + ft[0];
662     fy += fn[1] + ft[1];
663     fz += fn[2] + ft[2];
664   }
665 
666   f[i][0] += fx;
667   f[i][1] += fy;
668   f[i][2] += fz;
669 
670   if (newton_pair || j < nlocal) {
671     f[j][0] -= fx;
672     f[j][1] -= fy;
673     f[j][2] -= fz;
674   }
675 
676   if (evflag) ev_tally_xyz(i,j,nlocal,newton_pair,
677                            energy,0.0,fx,fy,fz,delx,dely,delz);
678 }
679 
680 /* ----------------------------------------------------------------------
681    Determine the interaction mode between i's vertices against j's edges
682 
683    i = atom i (body i)
684    j = atom j (body j)
685    x      = atoms' coordinates
686    f      = atoms' forces
687    torque = atoms' torques
688    tag    = atoms' tags
689    contact_list = list of contacts
690    num_contacts = number of contacts between i's vertices and j's edges
691    Return:
692      interact = 0 no interaction at all
693                 1 there's at least one case where i's vertices interacts
694                   with j's edges
695 ---------------------------------------------------------------------- */
696 
vertex_against_edge(int i,int j,double k_n,double k_na,double ** x,double ** f,double ** torque,tagint * tag,Contact * contact_list,int & num_contacts,double & evdwl,double * facc)697 int PairBodyRoundedPolygon::vertex_against_edge(int i, int j,
698                                                 double k_n, double k_na,
699                                                 double** x, double** f,
700                                                 double** torque, tagint* tag,
701                                                 Contact* contact_list,
702                                                 int &num_contacts,
703                                                 double &evdwl, double* facc)
704 {
705   int ni, npi, ifirst;
706   int nj, jfirst, nej, jefirst;
707   double xpi[3], xpj[3], dist, eradj, rradi, rradj;
708   double fx, fy, fz, energy;
709   int interact;
710 
711   npi = dnum[i];
712   ifirst = dfirst[i];
713   rradi = rounded_radius[i];
714 
715   jfirst = dfirst[j];
716   nej = ednum[j];
717   jefirst = edfirst[j];
718   eradj = enclosing_radius[j];
719   rradj = rounded_radius[j];
720 
721   energy = 0;
722   interact = 0;
723 
724   // loop through body i's vertices
725 
726   for (ni = 0; ni < npi; ni++) {
727 
728     // convert body-fixed coordinates to space-fixed, xi
729 
730     xpi[0] = x[i][0] + discrete[ifirst+ni][0];
731     xpi[1] = x[i][1] + discrete[ifirst+ni][1];
732     xpi[2] = x[i][2] + discrete[ifirst+ni][2];
733 
734     // compute the distance from the vertex to the COM of body j
735 
736     distance(xpi, x[j], dist);
737 
738     #ifdef _POLYGON_DEBUG
739     printf("Distance between vertex %d of body %d (%0.1f %0.1f %0.1f) "
740            "to body %d's COM: %f (cut = %0.1f)\n",
741            ni, xpi[0], xpi[1], xpi[2], atom->tag[i], atom->tag[j], dist,
742            eradj + rradi + rradj + cut_inner);
743     #endif
744 
745     // the vertex is within the enclosing circle (sphere) of body j,
746     // possibly interacting
747 
748     if (dist > eradj + rradj + rradi + cut_inner) continue;
749 
750     int mode, contact, p2vertex;
751     double d, R, hi[3], t, delx, dely, delz, fpair, shift;
752     double rij;
753 
754     // loop through body j's edges
755 
756     for (nj = 0; nj < nej; nj++) {
757 
758       // compute the distance between the edge nj to the vertex xpi
759 
760       mode = compute_distance_to_vertex(j, nj, x[j], rradj,
761                                         xpi, rradi, cut_inner,
762                                         d, hi, t, contact);
763 
764       if (mode == INVALID || mode == NONE) continue;
765 
766       if (mode == VERTEXI || mode == VERTEXJ) {
767 
768         interact = 1;
769 
770         // vertex i interacts with a vertex of the edge, but does not contact
771 
772         if (mode == VERTEXI) p2vertex = edge[jefirst+nj][0];
773         else if (mode == VERTEXJ) p2vertex = edge[jefirst+nj][1];
774 
775         // double xj[3];
776         // p2.body2space(p2vertex, xj);
777         xpj[0] = x[j][0] + discrete[jfirst+p2vertex][0];
778         xpj[1] = x[j][1] + discrete[jfirst+p2vertex][1];
779         xpj[2] = x[j][2] + discrete[jfirst+p2vertex][2];
780 
781         delx = xpi[0] - xpj[0];
782         dely = xpi[1] - xpj[1];
783         delz = xpi[2] - xpj[2];
784 
785         // R = surface separation = rij shifted by the rounded radii
786         // R = rij - (p1.rounded_radius + p2.rounded_radius);
787         // note: the force is defined for R, not for rij
788         // R > rc:     no interaction between vertex ni and p2vertex
789         // 0 < R < rc: cohesion between vertex ni and p2vertex
790         // R < 0:      deformation between vertex ni and p2vertex
791 
792         rij = sqrt(delx*delx + dely*dely + delz*delz);
793         R = rij - (rradi + rradj);
794         shift = k_na * cut_inner;
795 
796         // the normal frictional term -c_n * vn will be added later
797 
798         if (R <= 0) {           // deformation occurs
799           fpair = -k_n * R - shift;
800           energy += (0.5 * k_n * R + shift) * R;
801         } else if (R <= cut_inner) {   // not deforming but cohesive ranges overlap
802           fpair = k_na * R - shift;
803           energy += (-0.5 * k_na * R + shift) * R;
804         } else fpair = 0.0;
805 
806         fx = delx*fpair/rij;
807         fy = dely*fpair/rij;
808         fz = delz*fpair/rij;
809 
810         #ifdef _POLYGON_DEBUG
811         printf("  Interaction between vertex %d of %d and vertex %d of %d:",
812                ni, tag[i], p2vertex, tag[j]);
813         printf("    mode = %d; contact = %d; d = %f; rij = %f, t = %f\n",
814                mode, contact, d, rij, t);
815         printf("    R = %f; cut_inner = %f\n", R, cut_inner);
816         printf("    fpair = %f\n", fpair);
817         #endif
818 
819         // add forces to body i and body j directly
820         // avoid double counts this pair of vertices
821         // i and j can be either local or ghost atoms (bodies)
822         // probably need more work here when the vertices' interaction
823         // are not symmetric, e.g. j interacts with the edge
824         // consisting of i but in mode = EDGE instead of VERTEX*.
825         // OR, for the time being assume that the edge length is
826         // sufficiently greater than the rounded radius to distinguish
827         // vertex-vertex from vertex-edge contact modes.
828         // Special case: when i is a sphere, also accumulate
829 
830         if (tag[i] < tag[j] || npi == 1) {
831 
832           f[i][0] += fx;
833           f[i][1] += fy;
834           f[i][2] += fz;
835           sum_torque(x[i], xpi, fx, fy, fz, torque[i]);
836 
837           f[j][0] -= fx;
838           f[j][1] -= fy;
839           f[j][2] -= fz;
840           sum_torque(x[j], xpj, -fx, -fy, -fz, torque[j]);
841 
842           facc[0] += fx; facc[1] += fy; facc[2] += fz;
843 
844           #ifdef _POLYGON_DEBUG
845           printf("    from vertex-vertex: "
846                  "force on vertex %d of body %d: fx %f fy %f fz %f\n"
847                  "      torque body %d: %f %f %f\n"
848                  "      torque body %d: %f %f %f\n", ni, tag[i], fx, fy, fz,
849             tag[i],torque[i][0],torque[i][1],torque[i][2],
850             tag[j],torque[j][0],torque[j][1],torque[j][2]);
851         #endif
852         }
853 
854         #ifdef _CONVEX_POLYGON
855         // done with the edges from body j,
856         // given that vertex ni interacts with only one vertex
857         //   from one edge of body j
858         break;
859         #endif
860 
861       } else if (mode == EDGE) {
862 
863         interact = 1;
864 
865         // vertex i interacts with the edge
866 
867         delx = xpi[0] - hi[0];
868         dely = xpi[1] - hi[1];
869         delz = xpi[2] - hi[2];
870 
871         // R = surface separation = d shifted by the rounded radii
872         // R = d - (p1.rounded_radius + p2.rounded_radius);
873         // Note: the force is defined for R, not for d
874         // R > rc:     no interaction between vertex i and edge j
875         // 0 < R < rc: cohesion between vertex i and edge j
876         // R < 0:      deformation between vertex i and edge j
877         // rij = sqrt(delx*delx + dely*dely + delz*delz);
878 
879         R = d - (rradi + rradj);
880         shift = k_na * cut_inner;
881 
882         // the normal frictional term -c_n * vn will be added later
883 
884         if (R <= 0) {           // deformation occurs
885           fpair = -k_n * R - shift;
886           energy += (0.5 * k_n * R + shift) * R;
887         } else if (R <= cut_inner) {   // not deforming but cohesive ranges overlap
888           fpair = k_na * R - shift;
889           energy += (-0.5 * k_na * R + shift) * R;
890         } else fpair = 0.0;
891 
892         fx = delx*fpair/d;
893         fy = dely*fpair/d;
894         fz = delz*fpair/d;
895 
896         #ifdef _POLYGON_DEBUG
897         printf("  Interaction between vertex %d of %d and edge %d of %d:",
898                ni, tag[i], nj, tag[j]);
899         printf("    mode = %d; contact = %d; d = %f; t = %f\n",
900                mode, contact, d, t);
901         printf("    R = %f; cut_inner = %f\n", R, cut_inner);
902         printf("    fpair = %f\n", fpair);
903         #endif
904 
905         if (contact == 1) {
906 
907           // vertex ni of body i contacts with edge nj of body j
908 
909           contact_list[num_contacts].ibody = i;
910           contact_list[num_contacts].jbody = j;
911           contact_list[num_contacts].vertex = ni;
912           contact_list[num_contacts].edge = nj;
913           contact_list[num_contacts].xv[0] = xpi[0];
914           contact_list[num_contacts].xv[1] = xpi[1];
915           contact_list[num_contacts].xv[2] = xpi[2];
916           contact_list[num_contacts].xe[0] = hi[0];
917           contact_list[num_contacts].xe[1] = hi[1];
918           contact_list[num_contacts].xe[2] = hi[2];
919           contact_list[num_contacts].separation = R;
920           num_contacts++;
921 
922           // store forces to vertex ni and the edge nj
923           // to be rescaled later
924 
925           discrete[ifirst+ni][3] = fx;
926           discrete[ifirst+ni][4] = fy;
927           discrete[ifirst+ni][5] = fz;
928 
929           edge[jefirst+nj][2] = -fx;
930           edge[jefirst+nj][3] = -fy;
931           edge[jefirst+nj][4] = -fz;
932 
933           #ifdef _POLYGON_DEBUG
934           printf("  Stored forces at vertex and edge for accumulating later.\n");
935           #endif
936 
937         } else { // no contact
938 
939           // accumulate force and torque to both bodies directly
940 
941           f[i][0] += fx;
942           f[i][1] += fy;
943           f[i][2] += fz;
944           sum_torque(x[i], xpi, fx, fy, fz, torque[i]);
945 
946           f[j][0] -= fx;
947           f[j][1] -= fy;
948           f[j][2] -= fz;
949           sum_torque(x[j], hi, -fx, -fy, -fz, torque[j]);
950 
951           facc[0] += fx; facc[1] += fy; facc[2] += fz;
952 
953           #ifdef _POLYGON_DEBUG
954           printf("    from vertex-edge, no contact: "
955                  "force on vertex %d of body %d: fx %f fy %f fz %f\n"
956                  "      torque body %d: %f %f %f\n"
957                  "      torque body %d: %f %f %f\n", ni, tag[i], fx, fy, fz,
958                  tag[i],torque[i][0],torque[i][1],torque[i][2],
959                  tag[j],torque[j][0],torque[j][1],torque[j][2]);
960           #endif
961         } // end if contact
962 
963         #ifdef _CONVEX_POLYGON
964         // done with the edges from body j,
965         // given that vertex ni interacts with only one edge from body j
966         break;
967         #endif
968       } // end if mode
969 
970     } // end for looping through the edges of body j
971 
972   } // end for looping through the vertices of body i
973 
974   evdwl += energy;
975 
976   return interact;
977 }
978 
979 /* -------------------------------------------------------------------------
980   Compute the distance between an edge of body i and a vertex from
981   another body
982   Input:
983     ibody      = body i (i.e. atom i)
984     edge_index = edge index of body i
985     xmi        = atom i's coordinates (body i's center of mass)
986     x0         = coordinate of the tested vertex from another body
987     x0_rounded_radius = rounded radius of the tested vertex
988     cut_inner  = cutoff for vertex-vertex and vertex-edge interaction
989   Output:
990     d          = Distance from a point x0 to an edge
991     hi         = coordinates of the projection of x0 on the edge
992     t          = ratio to determine the relative position of hi
993                  wrt xi and xj on the segment
994   contact      = 0 no contact between the queried vertex and the edge
995                  1 contact detected
996   return
997     INVALID if the edge index is invalid
998     NONE    if there is no interaction
999     VERTEXI if the tested vertex interacts with the first vertex of the edge
1000     VERTEXJ if the tested vertex interacts with the second vertex of the edge
1001     EDGE    if the tested vertex interacts with the edge
1002 ------------------------------------------------------------------------- */
1003 
compute_distance_to_vertex(int ibody,int edge_index,double * xmi,double rounded_radius,double * x0,double x0_rounded_radius,double cut_inner,double & d,double hi[3],double & t,int & contact)1004 int PairBodyRoundedPolygon::compute_distance_to_vertex(int ibody,
1005                                                 int edge_index,
1006                                                 double *xmi,
1007                                                 double rounded_radius,
1008                                                 double* x0,
1009                                                 double x0_rounded_radius,
1010                                                 double cut_inner,
1011                                                 double &d,
1012                                                 double hi[3],
1013                                                 double &t,
1014                                                 int &contact)
1015 {
1016   if (edge_index >= ednum[ibody]) return INVALID;
1017 
1018   int mode,ifirst,iefirst,npi1,npi2;
1019   double xi1[3],xi2[3],u[3],v[3],uij[3];
1020   double udotv, magv, magucostheta;
1021   double delx,dely,delz;
1022 
1023   ifirst = dfirst[ibody];
1024   iefirst = edfirst[ibody];
1025   npi1 = static_cast<int>(edge[iefirst+edge_index][0]);
1026   npi2 = static_cast<int>(edge[iefirst+edge_index][1]);
1027 
1028   // compute the space-fixed coordinates for the vertices of the edge
1029 
1030   xi1[0] = xmi[0] + discrete[ifirst+npi1][0];
1031   xi1[1] = xmi[1] + discrete[ifirst+npi1][1];
1032   xi1[2] = xmi[2] + discrete[ifirst+npi1][2];
1033 
1034   xi2[0] = xmi[0] + discrete[ifirst+npi2][0];
1035   xi2[1] = xmi[1] + discrete[ifirst+npi2][1];
1036   xi2[2] = xmi[2] + discrete[ifirst+npi2][2];
1037 
1038   // u = x0 - xi1
1039 
1040   u[0] = x0[0] - xi1[0];
1041   u[1] = x0[1] - xi1[1];
1042   u[2] = x0[2] - xi1[2];
1043 
1044   // v = xi2 - xi1
1045 
1046   v[0] = xi2[0] - xi1[0];
1047   v[1] = xi2[1] - xi1[1];
1048   v[2] = xi2[2] - xi1[2];
1049 
1050   // dot product between u and v = magu * magv * costheta
1051 
1052   udotv = u[0] * v[0] + u[1] * v[1] + u[2] * v[2];
1053   magv = sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
1054   magucostheta = udotv / magv;
1055 
1056   // uij is the unit vector pointing from xi to xj
1057 
1058   uij[0] = v[0] / magv;
1059   uij[1] = v[1] / magv;
1060   uij[2] = v[2] / magv;
1061 
1062   // position of the projection of x0 on the line (xi, xj)
1063 
1064   hi[0] = xi1[0] + magucostheta * uij[0];
1065   hi[1] = xi1[1] + magucostheta * uij[1];
1066   hi[2] = xi1[2] + magucostheta * uij[2];
1067 
1068   // distance from x0 to the line (xi, xj) = distance from x0 to hi
1069 
1070   distance(hi, x0, d);
1071 
1072   // determine the interaction mode
1073   // for 2D: a vertex can interact with one edge at most
1074   // for 3D: a vertex can interact with one face at most
1075 
1076   mode = NONE;
1077   contact = 0;
1078 
1079   if (d > rounded_radius + x0_rounded_radius + cut_inner) {
1080 
1081     // if the vertex is far away from the edge
1082 
1083     mode = NONE;
1084 
1085   } else {
1086 
1087     // check if x0 (the queried vertex) and xmi (the body's center of mass)
1088     // are on the different sides of the edge
1089 
1090     #ifdef _CONVEX_POLYGON
1091     int m = opposite_sides(xi1, xi2, x0, xmi);
1092     #else
1093     int m = 1;
1094     #endif
1095 
1096     if (m == 0) {
1097 
1098       // x0 and xmi are on not the opposite sides of the edge
1099       // leave xpi for another edge to detect
1100 
1101       mode = NONE;
1102 
1103     } else {
1104 
1105       // x0 and xmi are on the different sides
1106       // t is the ratio to detect if x0 is closer to the vertices xi or xj
1107 
1108       if (fabs(xi2[0] - xi1[0]) > EPSILON)
1109         t = (hi[0] - xi1[0]) / (xi2[0] - xi1[0]);
1110       else if (fabs(xi2[1] - xi1[1]) > EPSILON)
1111         t = (hi[1] - xi1[1]) / (xi2[1] - xi1[1]);
1112       else if (fabs(xi2[2] - xi1[2]) > EPSILON)
1113         t = (hi[2] - xi1[2]) / (xi2[2] - xi1[2]);
1114 
1115       double contact_dist = rounded_radius + x0_rounded_radius;
1116       if (t >= 0 && t <= 1) {
1117         mode = EDGE;
1118         if (d < contact_dist + EPSILON)
1119           contact = 1;
1120 
1121       } else { // t < 0 || t > 1: closer to either vertices of the edge
1122 
1123         if (t < 0) {
1124           // measure the distance from x0 to xi1
1125           delx = x0[0] - xi1[0];
1126           dely = x0[1] - xi1[1];
1127           delz = x0[2] - xi1[2];
1128           double dx0xi1 = sqrt(delx*delx + dely*dely + delz*delz);
1129           if (dx0xi1 > contact_dist + cut_inner)
1130             mode = NONE;
1131           else
1132             mode = VERTEXI;
1133         } else {
1134           // measure the distance from x0 to xi2
1135           delx = x0[0] - xi2[0];
1136           dely = x0[1] - xi2[1];
1137           delz = x0[2] - xi2[2];
1138           double dx0xi2 = sqrt(delx*delx + dely*dely + delz*delz);
1139           if (dx0xi2 > contact_dist + cut_inner)
1140             mode = NONE;
1141           else
1142             mode = VERTEXJ;
1143         }
1144       } // end if t >= 0 && t <= 1
1145     } // end if x0 and xmi are on the same side of the edge
1146   }
1147 
1148   return mode;
1149 }
1150 
1151 /* ----------------------------------------------------------------------
1152   Compute contact forces between two bodies
1153   modify the force stored at the vertex and edge in contact by j_a
1154   sum forces and torque to the corresponding bodies
1155   fn = normal friction component
1156   ft = tangential friction component (-c_t * v_t)
1157 ------------------------------------------------------------------------- */
1158 
contact_forces(Contact & contact,double j_a,double ** x,double ** v,double ** angmom,double ** f,double ** torque,double &,double * facc)1159 void PairBodyRoundedPolygon::contact_forces(Contact& contact, double j_a,
1160                        double** x, double** v, double** angmom, double** f,
1161                        double** torque, double &/*evdwl*/, double* facc)
1162 {
1163   int ibody,jbody,ibonus,jbonus,ifirst,jefirst,ni,nj;
1164   double fx,fy,fz,delx,dely,delz,rsq,rsqinv;
1165   double vr1,vr2,vr3,vnnr,vn1,vn2,vn3,vt1,vt2,vt3;
1166   double fn[3],ft[3],vi[3],vj[3];
1167   double *quat, *inertia;
1168   AtomVecBody::Bonus *bonus;
1169 
1170   ibody = contact.ibody;
1171   jbody = contact.jbody;
1172 
1173   // compute the velocity of the vertex in the space-fixed frame
1174 
1175   ibonus = atom->body[ibody];
1176   bonus = &avec->bonus[ibonus];
1177   quat = bonus->quat;
1178   inertia = bonus->inertia;
1179   total_velocity(contact.xv, x[ibody], v[ibody], angmom[ibody],
1180                  inertia, quat, vi);
1181 
1182   // compute the velocity of the point on the edge in the space-fixed frame
1183 
1184   jbonus = atom->body[jbody];
1185   bonus = &avec->bonus[jbonus];
1186   quat = bonus->quat;
1187   inertia = bonus->inertia;
1188   total_velocity(contact.xe, x[jbody], v[jbody], angmom[jbody],
1189                  inertia, quat, vj);
1190 
1191   // vector pointing from the vertex to the point on the edge
1192 
1193   delx = contact.xv[0] - contact.xe[0];
1194   dely = contact.xv[1] - contact.xe[1];
1195   delz = contact.xv[2] - contact.xe[2];
1196   rsq = delx*delx + dely*dely + delz*delz;
1197   rsqinv = 1.0/rsq;
1198 
1199   // relative translational velocity
1200 
1201   vr1 = vi[0] - vj[0];
1202   vr2 = vi[1] - vj[1];
1203   vr3 = vi[2] - vj[2];
1204 
1205   // normal component
1206 
1207   vnnr = vr1*delx + vr2*dely + vr3*delz;
1208   vn1 = delx*vnnr * rsqinv;
1209   vn2 = dely*vnnr * rsqinv;
1210   vn3 = delz*vnnr * rsqinv;
1211 
1212   // tangential component
1213 
1214   vt1 = vr1 - vn1;
1215   vt2 = vr2 - vn2;
1216   vt3 = vr3 - vn3;
1217 
1218   // normal friction term at contact
1219 
1220   fn[0] = -c_n * vn1;
1221   fn[1] = -c_n * vn2;
1222   fn[2] = -c_n * vn3;
1223 
1224   // tangential friction term at contact
1225   // excluding the tangential deformation term for now
1226 
1227   ft[0] = -c_t * vt1;
1228   ft[1] = -c_t * vt2;
1229   ft[2] = -c_t * vt3;
1230 
1231   // only the cohesive force is scaled by j_a
1232   // mu * fne = tangential friction deformation during gross sliding
1233   // see Eq. 4, Fraige et al.
1234 
1235   ifirst = dfirst[ibody];
1236   ni = contact.vertex;
1237 
1238   fx = discrete[ifirst+ni][3] * j_a + fn[0] + ft[0] +
1239     mu * discrete[ifirst+ni][3];
1240   fy = discrete[ifirst+ni][4] * j_a + fn[1] + ft[1] +
1241     mu * discrete[ifirst+ni][4];
1242   fz = discrete[ifirst+ni][5] * j_a + fn[2] + ft[2] +
1243     mu * discrete[ifirst+ni][5];
1244   f[ibody][0] += fx;
1245   f[ibody][1] += fy;
1246   f[ibody][2] += fz;
1247   sum_torque(x[ibody], contact.xv, fx, fy, fz, torque[ibody]);
1248 
1249   // accumulate forces to the vertex only
1250 
1251   facc[0] += fx; facc[1] += fy; facc[2] += fz;
1252 
1253   // only the cohesive force is scaled by j_a
1254   // mu * fne = tangential friction deformation during gross sliding
1255   // Eq. 4, Fraige et al.
1256 
1257   jefirst = edfirst[jbody];
1258   nj = contact.edge;
1259 
1260   fx = edge[jefirst+nj][2] * j_a - fn[0] - ft[0] +
1261     mu * edge[jefirst+nj][2];
1262   fy = edge[jefirst+nj][3] * j_a - fn[1] - ft[1] +
1263     mu * edge[jefirst+nj][3];
1264   fz = edge[jefirst+nj][4] * j_a - fn[2] - ft[2] +
1265     mu * edge[jefirst+nj][4];
1266   f[jbody][0] += fx;
1267   f[jbody][1] += fy;
1268   f[jbody][2] += fz;
1269   sum_torque(x[jbody], contact.xe, fx, fy, fz, torque[jbody]);
1270 
1271   #ifdef _POLYGON_DEBUG
1272   printf("From contact forces: vertex fx %f fy %f fz %f\n"
1273          "      torque body %d: %f %f %f\n"
1274          "      torque body %d: %f %f %f\n",
1275          discrete[ifirst+ni][3], discrete[ifirst+ni][4], discrete[ifirst+ni][5],
1276          atom->tag[ibody],torque[ibody][0],torque[ibody][1],torque[ibody][2],
1277          atom->tag[jbody],torque[jbody][0],torque[jbody][1],torque[jbody][2]);
1278   #endif
1279 }
1280 
1281 /* ----------------------------------------------------------------------
1282   Determine the length of the contact segment, i.e. the separation between
1283   2 contacts, should be extended for 3D models.
1284 ------------------------------------------------------------------------- */
1285 
contact_separation(const Contact & c1,const Contact & c2)1286 double PairBodyRoundedPolygon::contact_separation(const Contact& c1,
1287                                                   const Contact& c2)
1288 {
1289   double x1 = c1.xv[0];
1290   double y1 = c1.xv[1];
1291   double x2 = c1.xe[0];
1292   double y2 = c1.xe[1];
1293   double x3 = c2.xv[0];
1294   double y3 = c2.xv[1];
1295 
1296   double delta_a = 0.0;
1297   if (fabs(x2 - x1) > EPSILON) {
1298     double A = (y2 - y1) / (x2 - x1);
1299     delta_a = fabs(y1 - A * x1 - y3 + A * x3) / sqrt(1 + A * A);
1300   } else {
1301     delta_a = fabs(x1 - x3);
1302   }
1303 
1304   return delta_a;
1305 }
1306 
1307 /* ----------------------------------------------------------------------
1308   Accumulate torque to body from the force f=(fx,fy,fz) acting at point x
1309 ------------------------------------------------------------------------- */
1310 
sum_torque(double * xm,double * x,double fx,double fy,double fz,double * torque)1311 void PairBodyRoundedPolygon::sum_torque(double* xm, double *x, double fx,
1312                                         double fy, double fz, double* torque)
1313 {
1314   double rx = x[0] - xm[0];
1315   double ry = x[1] - xm[1];
1316   double rz = x[2] - xm[2];
1317   double tx = ry * fz - rz * fy;
1318   double ty = rz * fx - rx * fz;
1319   double tz = rx * fy - ry * fx;
1320   torque[0] += tx;
1321   torque[1] += ty;
1322   torque[2] += tz;
1323 }
1324 
1325 /* ----------------------------------------------------------------------
1326   Test if two points a and b are in opposite sides of the line that
1327   connects two points x1 and x2
1328 ------------------------------------------------------------------------- */
1329 
opposite_sides(double * x1,double * x2,double * a,double * b)1330 int PairBodyRoundedPolygon::opposite_sides(double* x1, double* x2,
1331                                            double* a, double* b)
1332 {
1333   double m_a = (x1[1] - x2[1])*(a[0] - x1[0]) + (x2[0] - x1[0])*(a[1] - x1[1]);
1334   double m_b = (x1[1] - x2[1])*(b[0] - x1[0]) + (x2[0] - x1[0])*(b[1] - x1[1]);
1335   // equal to zero when either a or b is inline with the line x1-x2
1336   if (m_a * m_b <= 0)
1337     return 1;
1338   else
1339     return 0;
1340 }
1341 
1342 /* ----------------------------------------------------------------------
1343   Calculate the total velocity of a point (vertex, a point on an edge):
1344     vi = vcm + omega ^ (p - xcm)
1345 ------------------------------------------------------------------------- */
1346 
total_velocity(double * p,double * xcm,double * vcm,double * angmom,double * inertia,double * quat,double * vi)1347 void PairBodyRoundedPolygon::total_velocity(double* p, double *xcm,
1348                               double* vcm, double *angmom, double *inertia,
1349                               double *quat, double* vi)
1350 {
1351   double r[3],omega[3],ex_space[3],ey_space[3],ez_space[3];
1352   r[0] = p[0] - xcm[0];
1353   r[1] = p[1] - xcm[1];
1354   r[2] = p[2] - xcm[2];
1355   MathExtra::q_to_exyz(quat,ex_space,ey_space,ez_space);
1356   MathExtra::angmom_to_omega(angmom,ex_space,ey_space,ez_space,
1357                              inertia,omega);
1358   vi[0] = omega[1]*r[2] - omega[2]*r[1] + vcm[0];
1359   vi[1] = omega[2]*r[0] - omega[0]*r[2] + vcm[1];
1360   vi[2] = omega[0]*r[1] - omega[1]*r[0] + vcm[2];
1361 }
1362 
1363 /* ---------------------------------------------------------------------- */
1364 
distance(const double * x2,const double * x1,double & r)1365 void PairBodyRoundedPolygon::distance(const double* x2, const double* x1,
1366                                       double& r)
1367 {
1368   r = sqrt((x2[0] - x1[0]) * (x2[0] - x1[0])
1369     + (x2[1] - x1[1]) * (x2[1] - x1[1])
1370     + (x2[2] - x1[2]) * (x2[2] - x1[2]));
1371 }
1372