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 #include "region.h"
16
17 #include "domain.h"
18 #include "error.h"
19 #include "input.h"
20 #include "lattice.h"
21 #include "math_extra.h"
22 #include "update.h"
23 #include "variable.h"
24
25 #include <cmath>
26 #include <cstring>
27
28 using namespace LAMMPS_NS;
29
30 /* ---------------------------------------------------------------------- */
31
Region(LAMMPS * lmp,int,char ** arg)32 Region::Region(LAMMPS *lmp, int /*narg*/, char **arg) :
33 Pointers(lmp),
34 id(nullptr), style(nullptr), contact(nullptr), list(nullptr),
35 xstr(nullptr), ystr(nullptr), zstr(nullptr), tstr(nullptr)
36 {
37 id = utils::strdup(arg[0]);
38 style = utils::strdup(arg[1]);
39
40 varshape = 0;
41 xstr = ystr = zstr = tstr = nullptr;
42 dx = dy = dz = 0.0;
43
44 size_restart = 5;
45 Region::reset_vel();
46 copymode = 0;
47 list = nullptr;
48 nregion = 1;
49 }
50
51 /* ---------------------------------------------------------------------- */
52
~Region()53 Region::~Region()
54 {
55 if (copymode) return;
56
57 delete [] id;
58 delete [] style;
59
60 delete [] xstr;
61 delete [] ystr;
62 delete [] zstr;
63 delete [] tstr;
64 }
65
66 /* ---------------------------------------------------------------------- */
67
init()68 void Region::init()
69 {
70 if (xstr) {
71 xvar = input->variable->find(xstr);
72 if (xvar < 0) error->all(FLERR,"Variable name for region does not exist");
73 if (!input->variable->equalstyle(xvar))
74 error->all(FLERR,"Variable for region is invalid style");
75 }
76 if (ystr) {
77 yvar = input->variable->find(ystr);
78 if (yvar < 0) error->all(FLERR,"Variable name for region does not exist");
79 if (!input->variable->equalstyle(yvar))
80 error->all(FLERR,"Variable for region is not equal style");
81 }
82 if (zstr) {
83 zvar = input->variable->find(zstr);
84 if (zvar < 0) error->all(FLERR,"Variable name for region does not exist");
85 if (!input->variable->equalstyle(zvar))
86 error->all(FLERR,"Variable for region is not equal style");
87 }
88 if (tstr) {
89 tvar = input->variable->find(tstr);
90 if (tvar < 0) error->all(FLERR,"Variable name for region does not exist");
91 if (!input->variable->equalstyle(tvar))
92 error->all(FLERR,"Variable for region is not equal style");
93 }
94 vel_timestep = -1;
95 }
96
97 /* ----------------------------------------------------------------------
98 return 1 if region is dynamic (moves/rotates) or has variable shape
99 else return 0 if static
100 ------------------------------------------------------------------------- */
101
dynamic_check()102 int Region::dynamic_check()
103 {
104 if (dynamic || varshape) return 1;
105 return 0;
106 }
107
108 /* ----------------------------------------------------------------------
109 called before looping over atoms with match() or surface()
110 this insures any variables used by region are invoked once per timestep
111 also insures variables are invoked by all procs even those w/out atoms
112 necessary if equal-style variable invokes global operation
113 with MPI_Allreduce, e.g. xcm() or count()
114 ------------------------------------------------------------------------- */
115
prematch()116 void Region::prematch()
117 {
118 if (varshape) shape_update();
119 if (dynamic) pretransform();
120 }
121
122 /* ----------------------------------------------------------------------
123 determine if point x,y,z is a match to region volume
124 XOR computes 0 if 2 args are the same, 1 if different
125 note that inside() returns 1 for points on surface of region
126 thus point on surface of exterior region will not match
127 if region has variable shape, invoke shape_update() once per timestep
128 if region is dynamic, apply inverse transform to x,y,z
129 unmove first, then unrotate, so don't have to change rotation point
130 caller is responsible for wrapping this call with
131 modify->clearstep_compute() and modify->addstep_compute() if needed
132 ------------------------------------------------------------------------- */
133
match(double x,double y,double z)134 int Region::match(double x, double y, double z)
135 {
136 if (dynamic) inverse_transform(x,y,z);
137 if (openflag) return 1;
138 return !(inside(x,y,z) ^ interior);
139 }
140
141 /* ----------------------------------------------------------------------
142 generate list of contact points for interior or exterior regions
143 if region has variable shape, invoke shape_update() once per timestep
144 if region is dynamic:
145 before: inverse transform x,y,z (unmove, then unrotate)
146 after: forward transform contact point xs,yx,zs (rotate, then move),
147 then reset contact delx,dely,delz based on new contact point
148 no need to do this if no rotation since delxyz doesn't change
149 caller is responsible for wrapping this call with
150 modify->clearstep_compute() and modify->addstep_compute() if needed
151 ------------------------------------------------------------------------- */
152
surface(double x,double y,double z,double cutoff)153 int Region::surface(double x, double y, double z, double cutoff)
154 {
155 int ncontact;
156 double xs,ys,zs;
157 double xnear[3],xorig[3];
158
159 if (dynamic) {
160 xorig[0] = x;
161 xorig[1] = y;
162 xorig[2] = z;
163 inverse_transform(x,y,z);
164 }
165
166 xnear[0] = x;
167 xnear[1] = y;
168 xnear[2] = z;
169
170 if (!openflag) {
171 if (interior) ncontact = surface_interior(xnear,cutoff);
172 else ncontact = surface_exterior(xnear,cutoff);
173 } else {
174 // one of surface_int/ext() will return 0
175 // so no need to worry about offset of contact indices
176 ncontact = surface_exterior(xnear,cutoff) + surface_interior(xnear,cutoff);
177 }
178
179 if (rotateflag && ncontact) {
180 for (int i = 0; i < ncontact; i++) {
181 xs = xnear[0] - contact[i].delx;
182 ys = xnear[1] - contact[i].dely;
183 zs = xnear[2] - contact[i].delz;
184 forward_transform(xs,ys,zs);
185 contact[i].delx = xorig[0] - xs;
186 contact[i].dely = xorig[1] - ys;
187 contact[i].delz = xorig[2] - zs;
188 }
189 }
190
191 return ncontact;
192 }
193
194 /* ----------------------------------------------------------------------
195 add a single contact at Nth location in contact array
196 x = particle position
197 xp,yp,zp = region surface point
198 ------------------------------------------------------------------------- */
199
add_contact(int n,double * x,double xp,double yp,double zp)200 void Region::add_contact(int n, double *x, double xp, double yp, double zp)
201 {
202 double delx = x[0] - xp;
203 double dely = x[1] - yp;
204 double delz = x[2] - zp;
205 contact[n].r = sqrt(delx*delx + dely*dely + delz*delz);
206 contact[n].radius = 0;
207 contact[n].delx = delx;
208 contact[n].dely = dely;
209 contact[n].delz = delz;
210 }
211
212 /* ----------------------------------------------------------------------
213 pre-compute dx,dy,dz and theta for a moving/rotating region
214 called once for the region before per-atom loop, via prematch()
215 ------------------------------------------------------------------------- */
216
pretransform()217 void Region::pretransform()
218 {
219 if (moveflag) {
220 if (xstr) dx = input->variable->compute_equal(xvar);
221 if (ystr) dy = input->variable->compute_equal(yvar);
222 if (zstr) dz = input->variable->compute_equal(zvar);
223 }
224 if (rotateflag) theta = input->variable->compute_equal(tvar);
225 }
226
227 /* ----------------------------------------------------------------------
228 transform a point x,y,z in region space to moved space
229 rotate first (around original P), then displace
230 ------------------------------------------------------------------------- */
231
forward_transform(double & x,double & y,double & z)232 void Region::forward_transform(double &x, double &y, double &z)
233 {
234 if (rotateflag) rotate(x,y,z,theta);
235 if (moveflag) {
236 x += dx;
237 y += dy;
238 z += dz;
239 }
240 }
241
242 /* ----------------------------------------------------------------------
243 transform a point x,y,z in moved space back to region space
244 undisplace first, then unrotate (around original P)
245 ------------------------------------------------------------------------- */
246
inverse_transform(double & x,double & y,double & z)247 void Region::inverse_transform(double &x, double &y, double &z)
248 {
249 if (moveflag) {
250 x -= dx;
251 y -= dy;
252 z -= dz;
253 }
254 if (rotateflag) rotate(x,y,z,-theta);
255 }
256
257 /* ----------------------------------------------------------------------
258 rotate x,y,z by angle via right-hand rule around point and runit normal
259 sign of angle determines whether rotating forward/backward in time
260 return updated x,y,z
261 R = vector axis of rotation
262 P = point = point to rotate around
263 R0 = runit = unit vector for R
264 X0 = x,y,z = initial coord of atom
265 D = X0 - P = vector from P to X0
266 C = (D dot R0) R0 = projection of D onto R, i.e. Dparallel
267 A = D - C = vector from R line to X0, i.e. Dperp
268 B = R0 cross A = vector perp to A in plane of rotation, same len as A
269 A,B define plane of circular rotation around R line
270 new x,y,z = P + C + A cos(angle) + B sin(angle)
271 ------------------------------------------------------------------------- */
272
rotate(double & x,double & y,double & z,double angle)273 void Region::rotate(double &x, double &y, double &z, double angle)
274 {
275 double a[3],b[3],c[3],d[3],disp[3];
276
277 double sine = sin(angle);
278 double cosine = cos(angle);
279 d[0] = x - point[0];
280 d[1] = y - point[1];
281 d[2] = z - point[2];
282 double x0dotr = d[0]*runit[0] + d[1]*runit[1] + d[2]*runit[2];
283 c[0] = x0dotr * runit[0];
284 c[1] = x0dotr * runit[1];
285 c[2] = x0dotr * runit[2];
286 a[0] = d[0] - c[0];
287 a[1] = d[1] - c[1];
288 a[2] = d[2] - c[2];
289 b[0] = runit[1]*a[2] - runit[2]*a[1];
290 b[1] = runit[2]*a[0] - runit[0]*a[2];
291 b[2] = runit[0]*a[1] - runit[1]*a[0];
292 disp[0] = a[0]*cosine + b[0]*sine;
293 disp[1] = a[1]*cosine + b[1]*sine;
294 disp[2] = a[2]*cosine + b[2]*sine;
295 x = point[0] + c[0] + disp[0];
296 y = point[1] + c[1] + disp[1];
297 z = point[2] + c[2] + disp[2];
298 }
299
300 /* ----------------------------------------------------------------------
301 parse optional parameters at end of region input line
302 ------------------------------------------------------------------------- */
303
options(int narg,char ** arg)304 void Region::options(int narg, char **arg)
305 {
306 if (narg < 0) error->all(FLERR,"Illegal region command");
307
308 // option defaults
309
310 interior = 1;
311 scaleflag = 1;
312 moveflag = rotateflag = 0;
313
314 openflag = 0;
315 for (int i = 0; i < 6; i++) open_faces[i] = 0;
316
317 int iarg = 0;
318 while (iarg < narg) {
319 if (strcmp(arg[iarg],"units") == 0) {
320 if (iarg+2 > narg) error->all(FLERR,"Illegal region command");
321 if (strcmp(arg[iarg+1],"box") == 0) scaleflag = 0;
322 else if (strcmp(arg[iarg+1],"lattice") == 0) scaleflag = 1;
323 else error->all(FLERR,"Illegal region command");
324 iarg += 2;
325 } else if (strcmp(arg[iarg],"side") == 0) {
326 if (iarg+2 > narg) error->all(FLERR,"Illegal region command");
327 if (strcmp(arg[iarg+1],"in") == 0) interior = 1;
328 else if (strcmp(arg[iarg+1],"out") == 0) interior = 0;
329 else error->all(FLERR,"Illegal region command");
330 iarg += 2;
331
332 } else if (strcmp(arg[iarg],"move") == 0) {
333 if (iarg+4 > narg) error->all(FLERR,"Illegal region command");
334 if (strcmp(arg[iarg+1],"NULL") != 0) {
335 if (strstr(arg[iarg+1],"v_") != arg[iarg+1])
336 error->all(FLERR,"Illegal region command");
337 xstr = utils::strdup(&arg[iarg+1][2]);
338 }
339 if (strcmp(arg[iarg+2],"NULL") != 0) {
340 if (strstr(arg[iarg+2],"v_") != arg[iarg+2])
341 error->all(FLERR,"Illegal region command");
342 ystr = utils::strdup(&arg[iarg+2][2]);
343 }
344 if (strcmp(arg[iarg+3],"NULL") != 0) {
345 if (strstr(arg[iarg+3],"v_") != arg[iarg+3])
346 error->all(FLERR,"Illegal region command");
347 zstr = utils::strdup(&arg[iarg+3][2]);
348 }
349 moveflag = 1;
350 iarg += 4;
351
352 } else if (strcmp(arg[iarg],"rotate") == 0) {
353 if (iarg+8 > narg) error->all(FLERR,"Illegal region command");
354 if (strstr(arg[iarg+1],"v_") != arg[iarg+1])
355 error->all(FLERR,"Illegal region command");
356 tstr = utils::strdup(&arg[iarg+1][2]);
357 point[0] = utils::numeric(FLERR,arg[iarg+2],false,lmp);
358 point[1] = utils::numeric(FLERR,arg[iarg+3],false,lmp);
359 point[2] = utils::numeric(FLERR,arg[iarg+4],false,lmp);
360 axis[0] = utils::numeric(FLERR,arg[iarg+5],false,lmp);
361 axis[1] = utils::numeric(FLERR,arg[iarg+6],false,lmp);
362 axis[2] = utils::numeric(FLERR,arg[iarg+7],false,lmp);
363 rotateflag = 1;
364 iarg += 8;
365
366 } else if (strcmp(arg[iarg],"open") == 0) {
367 if (iarg+2 > narg) error->all(FLERR,"Illegal region command");
368 int iface = utils::inumeric(FLERR,arg[iarg+1],false,lmp);
369 if (iface < 1 || iface > 6) error->all(FLERR,"Illegal region command");
370 // additional checks on valid face index are done by region classes
371 open_faces[iface-1] = 1;
372 openflag = 1;
373 iarg += 2;
374 }
375 else error->all(FLERR,"Illegal region command");
376 }
377
378 // error check
379
380 if ((moveflag || rotateflag) &&
381 (strcmp(style,"union") == 0 || strcmp(style,"intersect") == 0))
382 error->all(FLERR,"Region union or intersect cannot be dynamic");
383
384 // setup scaling
385
386 if (scaleflag) {
387 xscale = domain->lattice->xlattice;
388 yscale = domain->lattice->ylattice;
389 zscale = domain->lattice->zlattice;
390 }
391 else xscale = yscale = zscale = 1.0;
392
393 if (rotateflag) {
394 point[0] *= xscale;
395 point[1] *= yscale;
396 point[2] *= zscale;
397 }
398
399 // runit = unit vector along rotation axis
400
401 if (rotateflag) {
402 double len = sqrt(axis[0]*axis[0] + axis[1]*axis[1] + axis[2]*axis[2]);
403 if (len == 0.0)
404 error->all(FLERR,"Region cannot have 0 length rotation vector");
405 runit[0] = axis[0]/len;
406 runit[1] = axis[1]/len;
407 runit[2] = axis[2]/len;
408 }
409
410 if (moveflag || rotateflag) dynamic = 1;
411 else dynamic = 0;
412 }
413
414 /* ----------------------------------------------------------------------
415 find nearest point to C on line segment A,B and return it as D
416 project (C-A) onto (B-A)
417 t = length of that projection, normalized by length of (B-A)
418 t <= 0, C is closest to A
419 t >= 1, C is closest to B
420 else closest point is between A and B
421 ------------------------------------------------------------------------- */
422
point_on_line_segment(double * a,double * b,double * c,double * d)423 void Region::point_on_line_segment(double *a, double *b,
424 double *c, double *d)
425 {
426 double ba[3],ca[3];
427
428 MathExtra::sub3(b,a,ba);
429 MathExtra::sub3(c,a,ca);
430 double t = MathExtra::dot3(ca,ba) / MathExtra::dot3(ba,ba);
431 if (t <= 0.0) {
432 d[0] = a[0];
433 d[1] = a[1];
434 d[2] = a[2];
435 } else if (t >= 1.0) {
436 d[0] = b[0];
437 d[1] = b[1];
438 d[2] = b[2];
439 } else {
440 d[0] = a[0] + t*ba[0];
441 d[1] = a[1] + t*ba[1];
442 d[2] = a[2] + t*ba[2];
443 }
444 }
445
446 /* ----------------------------------------------------------------------
447 infer translational and angular velocity of region
448 necessary b/c motion variables are for displacement & theta
449 there is no analytic formula for v & omega
450 prev[4] contains values of dx,dy,dz,theta at previous step
451 used for difference, then updated to current step values
452 dt is time elapsed since previous step
453 rpoint = point updated by current displacement
454 called by fix wall/gran/region every timestep
455 ------------------------------------------------------------------------- */
456
set_velocity()457 void Region::set_velocity()
458 {
459 if (vel_timestep == update->ntimestep) return;
460 vel_timestep = update->ntimestep;
461 if (moveflag) {
462 if (update->ntimestep > 0) {
463 v[0] = (dx - prev[0])/update->dt;
464 v[1] = (dy - prev[1])/update->dt;
465 v[2] = (dz - prev[2])/update->dt;
466 }
467 else v[0] = v[1] = v[2] = 0.0;
468 prev[0] = dx;
469 prev[1] = dy;
470 prev[2] = dz;
471 }
472
473 if (rotateflag) {
474 rpoint[0] = point[0] + dx;
475 rpoint[1] = point[1] + dy;
476 rpoint[2] = point[2] + dz;
477 if (update->ntimestep > 0) {
478 double angvel = (theta-prev[3]) / update->dt;
479 omega[0] = angvel*axis[0];
480 omega[1] = angvel*axis[1];
481 omega[2] = angvel*axis[2];
482 }
483 else omega[0] = omega[1] = omega[2] = 0.0;
484 prev[3] = theta;
485 }
486
487 if (varshape) {
488 set_velocity_shape();
489 }
490 }
491
492 /* ----------------------------------------------------------------------
493 compute velocity of wall for given contact
494 since contacts only store delx/y/z, need to pass particle coords
495 to compute contact point
496 called by fix/wall/gran/region every contact every timestep
497 ------------------------------------------------------------------------- */
498
velocity_contact(double * vwall,double * x,int ic)499 void Region::velocity_contact(double *vwall, double *x, int ic)
500 {
501 double xc[3];
502
503 vwall[0] = vwall[1] = vwall[2] = 0.0;
504
505 if (moveflag) {
506 vwall[0] = v[0];
507 vwall[1] = v[1];
508 vwall[2] = v[2];
509 }
510 if (rotateflag) {
511 xc[0] = x[0] - contact[ic].delx;
512 xc[1] = x[1] - contact[ic].dely;
513 xc[2] = x[2] - contact[ic].delz;
514 vwall[0] += omega[1]*(xc[2] - rpoint[2]) - omega[2]*(xc[1] - rpoint[1]);
515 vwall[1] += omega[2]*(xc[0] - rpoint[0]) - omega[0]*(xc[2] - rpoint[2]);
516 vwall[2] += omega[0]*(xc[1] - rpoint[1]) - omega[1]*(xc[0] - rpoint[0]);
517 }
518
519 if (varshape && contact[ic].varflag) velocity_contact_shape(vwall, xc);
520 }
521
522
523 /* ----------------------------------------------------------------------
524 increment length of restart buffer based on region info
525 used by restart of fix/wall/gran/region
526 ------------------------------------------------------------------------- */
527
length_restart_string(int & n)528 void Region::length_restart_string(int &n)
529 {
530 n += sizeof(int) + strlen(id)+1 +
531 sizeof(int) + strlen(style)+1 + sizeof(int) +
532 size_restart*sizeof(double);
533 }
534
535 /* ----------------------------------------------------------------------
536 region writes its current style, id, number of sub-regions, position/angle
537 needed by fix/wall/gran/region to compute velocity by differencing scheme
538 ------------------------------------------------------------------------- */
539
write_restart(FILE * fp)540 void Region::write_restart(FILE *fp)
541 {
542 int sizeid = (strlen(id)+1);
543 int sizestyle = (strlen(style)+1);
544 fwrite(&sizeid, sizeof(int), 1, fp);
545 fwrite(id,1,sizeid,fp);
546 fwrite(&sizestyle,sizeof(int),1,fp);
547 fwrite(style,1,sizestyle,fp);
548 fwrite(&nregion,sizeof(int),1,fp);
549 fwrite(prev,sizeof(double),size_restart,fp);
550 }
551
552 /* ----------------------------------------------------------------------
553 region reads style, id, number of sub-regions from restart file
554 if they match current region, also read previous position/angle
555 needed by fix/wall/gran/region to compute velocity by differencing scheme
556 ------------------------------------------------------------------------- */
557
restart(char * buf,int & n)558 int Region::restart(char *buf, int &n)
559 {
560 int size = *((int *) (&buf[n]));
561 n += sizeof(int);
562 if ((size <= 0) || (strcmp(&buf[n],id) != 0)) return 0;
563 n += size;
564
565 size = *((int *) (&buf[n]));
566 n += sizeof(int);
567 if ((size <= 0) || (strcmp(&buf[n],style) != 0)) return 0;
568 n += size;
569
570 int restart_nreg = *((int *) (&buf[n]));
571 n += sizeof(int);
572 if (restart_nreg != nregion) return 0;
573
574 memcpy(prev,&buf[n],size_restart*sizeof(double));
575 return 1;
576 }
577
578 /* ----------------------------------------------------------------------
579 set prev vector to zero
580 ------------------------------------------------------------------------- */
581
reset_vel()582 void Region::reset_vel()
583 {
584 for (int i = 0; i < size_restart; i++) prev[i] = 0;
585 }
586