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