1 /* ----------------------------------------------------------------------
2 This is the
3
4 ██╗ ██╗ ██████╗ ██████╗ ██████╗ ██╗ ██╗████████╗███████╗
5 ██║ ██║██╔════╝ ██╔════╝ ██╔════╝ ██║ ██║╚══██╔══╝██╔════╝
6 ██║ ██║██║ ███╗██║ ███╗██║ ███╗███████║ ██║ ███████╗
7 ██║ ██║██║ ██║██║ ██║██║ ██║██╔══██║ ██║ ╚════██║
8 ███████╗██║╚██████╔╝╚██████╔╝╚██████╔╝██║ ██║ ██║ ███████║
9 ╚══════╝╚═╝ ╚═════╝ ╚═════╝ ╚═════╝ ╚═╝ ╚═╝ ╚═╝ ╚══════╝®
10
11 DEM simulation engine, released by
12 DCS Computing Gmbh, Linz, Austria
13 http://www.dcs-computing.com, office@dcs-computing.com
14
15 LIGGGHTS® is part of CFDEM®project:
16 http://www.liggghts.com | http://www.cfdem.com
17
18 Core developer and main author:
19 Christoph Kloss, christoph.kloss@dcs-computing.com
20
21 LIGGGHTS® is open-source, distributed under the terms of the GNU Public
22 License, version 2 or later. It is distributed in the hope that it will
23 be useful, but WITHOUT ANY WARRANTY; without even the implied warranty
24 of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. You should have
25 received a copy of the GNU General Public License along with LIGGGHTS®.
26 If not, see http://www.gnu.org/licenses . See also top-level README
27 and LICENSE files.
28
29 LIGGGHTS® and CFDEM® are registered trade marks of DCS Computing GmbH,
30 the producer of the LIGGGHTS® software and the CFDEM®coupling software
31 See http://www.cfdem.com/terms-trademark-policy for details.
32
33 -------------------------------------------------------------------------
34 Contributing author and copyright for this file:
35 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include <cmath>
47 #include <stdlib.h>
48 #include <string.h>
49 #include "region_sphere.h"
50 #include "update.h"
51 #include "input.h"
52 #include "variable.h"
53 #include "error.h"
54 #include "force.h"
55
56 using namespace LAMMPS_NS;
57
58 enum{CONSTANT,VARIABLE};
59
60 /* ---------------------------------------------------------------------- */
61
RegSphere(LAMMPS * lmp,int narg,char ** arg)62 RegSphere::RegSphere(LAMMPS *lmp, int narg, char **arg) :
63 Region(lmp, narg, arg)
64 {
65 options(narg-6,&arg[6]);
66
67 xc = xscale*force->numeric(FLERR,arg[2]);
68 yc = yscale*force->numeric(FLERR,arg[3]);
69 zc = zscale*force->numeric(FLERR,arg[4]);
70
71 rstr = NULL;
72 if (strstr(arg[5],"v_") == arg[5]) {
73 int n = strlen(&arg[5][2]) + 1;
74 rstr = new char[n];
75 strcpy(rstr,&arg[5][2]);
76 radius = 0.0;
77 rstyle = VARIABLE;
78 varshape = 1;
79 variable_check();
80 shape_update();
81 } else {
82 radius = xscale*force->numeric(FLERR,arg[5]);
83 rstyle = CONSTANT;
84 }
85
86 // error check
87
88 if (radius < 0.0) error->all(FLERR,"Illegal region sphere command");
89
90 // extent of sphere
91 // for variable radius, uses initial radius
92
93 if (interior) {
94 bboxflag = 1;
95 extent_xlo = xc - radius;
96 extent_xhi = xc + radius;
97 extent_ylo = yc - radius;
98 extent_yhi = yc + radius;
99 extent_zlo = zc - radius;
100 extent_zhi = zc + radius;
101 } else bboxflag = 0;
102
103 cmax = 1;
104 contact = new Contact[cmax];
105 }
106
107 /* ---------------------------------------------------------------------- */
108
~RegSphere()109 RegSphere::~RegSphere()
110 {
111 delete [] rstr;
112 delete [] contact;
113 }
114
115 /* ---------------------------------------------------------------------- */
116
init()117 void RegSphere::init()
118 {
119 Region::init();
120 if (rstr) variable_check();
121 }
122
123 /* ----------------------------------------------------------------------
124 inside = 1 if x,y,z is inside or on surface
125 inside = 0 if x,y,z is outside and not on surface
126 ------------------------------------------------------------------------- */
127
inside(double x,double y,double z)128 int RegSphere::inside(double x, double y, double z)
129 {
130 double delx = x - xc;
131 double dely = y - yc;
132 double delz = z - zc;
133 double r = sqrt(delx*delx + dely*dely + delz*delz);
134
135 if (r <= radius) return 1;
136 return 0;
137 }
138
139 /* ----------------------------------------------------------------------
140 one contact if 0 <= x < cutoff from inner surface of sphere
141 no contact if outside (possible if called from union/intersect)
142 delxyz = vector from nearest point on sphere to x
143 special case: no contact if x is at center of sphere
144 ------------------------------------------------------------------------- */
145
surface_interior(double * x,double cutoff)146 int RegSphere::surface_interior(double *x, double cutoff)
147 {
148 double delx = x[0] - xc;
149 double dely = x[1] - yc;
150 double delz = x[2] - zc;
151 double r = sqrt(delx*delx + dely*dely + delz*delz);
152 if (r > radius || r == 0.0) return 0;
153
154 double delta = radius - r;
155 if (delta < cutoff) {
156 contact[0].r = delta;
157 contact[0].delx = delx*(1.0-radius/r);
158 contact[0].dely = dely*(1.0-radius/r);
159 contact[0].delz = delz*(1.0-radius/r);
160 return 1;
161 }
162 return 0;
163 }
164
165 /* ----------------------------------------------------------------------
166 one contact if 0 <= x < cutoff from outer surface of sphere
167 no contact if inside (possible if called from union/intersect)
168 delxyz = vector from nearest point on sphere to x
169 ------------------------------------------------------------------------- */
170
surface_exterior(double * x,double cutoff)171 int RegSphere::surface_exterior(double *x, double cutoff)
172 {
173 double delx = x[0] - xc;
174 double dely = x[1] - yc;
175 double delz = x[2] - zc;
176 double r = sqrt(delx*delx + dely*dely + delz*delz);
177 if (r < radius) return 0;
178
179 double delta = r - radius;
180 if (delta < cutoff) {
181 contact[0].r = delta;
182 contact[0].delx = delx*(1.0-radius/r);
183 contact[0].dely = dely*(1.0-radius/r);
184 contact[0].delz = delz*(1.0-radius/r);
185 return 1;
186 }
187 return 0;
188 }
189
190 /* ----------------------------------------------------------------------
191 change region shape via variable evaluation
192 ------------------------------------------------------------------------- */
193
shape_update()194 void RegSphere::shape_update()
195 {
196 radius = xscale * input->variable->compute_equal(rvar);
197 if (radius < 0.0)
198 error->one(FLERR,"Variable evaluation in region gave bad value");
199 }
200
201 /* ----------------------------------------------------------------------
202 error check on existence of variable
203 ------------------------------------------------------------------------- */
204
variable_check()205 void RegSphere::variable_check()
206 {
207 rvar = input->variable->find(rstr);
208 if (rvar < 0)
209 error->all(FLERR,"Variable name for region sphere does not exist");
210 if (!input->variable->equalstyle(rvar))
211 error->all(FLERR,"Variable for region sphere is invalid style");
212 }
213