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 Andreas Aigner (JKU Linz)
36
37 Copyright 2009-2012 JKU Linz
38 ------------------------------------------------------------------------- */
39
40 #include <cmath>
41 #include <stdlib.h>
42 #include <string.h>
43 #include "atom.h"
44 #include "force.h"
45 #include "pair.h"
46 #include "modify.h"
47 #include "memory.h"
48 #include "domain.h"
49 #include "respa.h"
50 #include "update.h"
51 #include "error.h"
52 #include "sph_kernels.h"
53 #include "fix_wall_sph.h"
54
55 // include last to ensure correct macros
56 #include "domain_definitions.h"
57
58 using namespace LAMMPS_NS;
59 using namespace FixConst;
60
61 enum{XPLANE,YPLANE,ZPLANE,ZCYLINDER}; // XYZ PLANE need to be 0,1,2
62
63 /* ---------------------------------------------------------------------- */
64
FixWallSph(LAMMPS * lmp,int narg,char ** arg)65 FixWallSph::FixWallSph(LAMMPS *lmp, int narg, char **arg) :
66 FixSph(lmp, narg, arg)
67 {
68 // wallstyle args
69
70 int iarg = 3;
71 if (strcmp(arg[iarg],"xplane") == 0) {
72 if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/sph command");
73 wallstyle = XPLANE;
74 if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
75 else lo = force->numeric(FLERR,arg[iarg+1]);
76 if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
77 else hi = force->numeric(FLERR,arg[iarg+2]);
78 iarg += 3;
79 } else if (strcmp(arg[iarg],"yplane") == 0) {
80 if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/sph command");
81 wallstyle = YPLANE;
82 if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
83 else lo = force->numeric(FLERR,arg[iarg+1]);
84 if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
85 else hi = force->numeric(FLERR,arg[iarg+2]);
86 iarg += 3;
87 } else if (strcmp(arg[iarg],"zplane") == 0) {
88 if (narg < iarg+3) error->all(FLERR,"Illegal fix wall/sph command");
89 wallstyle = ZPLANE;
90 if (strcmp(arg[iarg+1],"NULL") == 0) lo = -BIG;
91 else lo = force->numeric(FLERR,arg[iarg+1]);
92 if (strcmp(arg[iarg+2],"NULL") == 0) hi = BIG;
93 else hi = force->numeric(FLERR,arg[iarg+2]);
94 iarg += 3;
95 } else if (strcmp(arg[iarg],"zcylinder") == 0) {
96 if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/gran command");
97 wallstyle = ZCYLINDER;
98 lo = hi = 0.0;
99 cylradius = force->numeric(FLERR,arg[iarg+1]);
100 iarg += 2;
101 }
102
103 // parameters for penetration force
104 if (narg < iarg+2) error->all(FLERR,"Illegal fix wall/sph command, not enough arguments for penetration force");
105 r0 = force->numeric(FLERR,arg[iarg]);
106 D = force->numeric(FLERR,arg[iarg+1]);
107 iarg += 2;
108
109 if (wallstyle == XPLANE && domain->xperiodic)
110 error->all(FLERR,"Cannot use wall in periodic dimension");
111 if (wallstyle == YPLANE && domain->yperiodic)
112 error->all(FLERR,"Cannot use wall in periodic dimension");
113 if (wallstyle == ZPLANE && domain->zperiodic)
114 error->all(FLERR,"Cannot use wall in periodic dimension");
115 if (wallstyle == ZCYLINDER && (domain->xperiodic || domain->yperiodic))
116 error->all(FLERR,"Cannot use wall in periodic dimension");
117
118 }
119
120 /* ---------------------------------------------------------------------- */
121
~FixWallSph()122 FixWallSph::~FixWallSph()
123 {
124
125 }
126
127 /* ---------------------------------------------------------------------- */
128
setmask()129 int FixWallSph::setmask()
130 {
131 int mask = 0;
132 mask |= POST_FORCE;
133 mask |= POST_FORCE_RESPA;
134 return mask;
135 }
136
137 /* ---------------------------------------------------------------------- */
138
init()139 void FixWallSph::init()
140 {
141 FixSph::init();
142
143 if (strcmp(update->integrate_style,"respa") == 0)
144 nlevels_respa = ((Respa *) update->integrate)->nlevels;
145 }
146
147 /* ---------------------------------------------------------------------- */
148
setup(int vflag)149 void FixWallSph::setup(int vflag)
150 {
151 if (strstr(update->integrate_style,"verlet"))
152 post_force(vflag);
153 else {
154 ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa-1);
155 post_force_respa(vflag,nlevels_respa-1,0);
156 ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa-1);
157 }
158 }
159
160 /* ---------------------------------------------------------------------- */
161
post_force(int vflag)162 void FixWallSph::post_force(int vflag)
163 {
164 double dx,dy,dz,del1,del2,delxy,delr,rsq,r,rinv;
165 double fwall;
166
167 double wlo = lo;
168 double whi = hi;
169
170 double frac,frac2; // for penetration force
171
172 // loop over all my atoms
173 // rsq = distance from wall
174 // dx,dy,dz = signed distance from wall
175 // skip atom if not close enough to wall
176 // if wall was set to NULL, it's skipped since lo/hi are infinity
177 // compute force on atom if close enough to wall
178 // via wall potential matched to pair potential
179
180 double **x = atom->x;
181 double **f = atom->f;
182 int *mask = atom->mask;
183 int nlocal = atom->nlocal;
184
185 for (int i = 0; i < nlocal; i++) {
186 if (mask[i] & groupbit) {
187
188 dx = dy = dz = 0.0;
189
190 if (wallstyle == XPLANE) {
191 del1 = x[i][0] - wlo;
192 del2 = whi - x[i][0];
193 if (del1 < del2) dx = del1;
194 else dx = -del2;
195 } else if (wallstyle == YPLANE) {
196 del1 = x[i][1] - wlo;
197 del2 = whi - x[i][1];
198 if (del1 < del2) dy = del1;
199 else dy = -del2;
200 } else if (wallstyle == ZPLANE) {
201 del1 = x[i][2] - wlo;
202 del2 = whi - x[i][2];
203 if (del1 < del2) dz = del1;
204 else dz = -del2;
205 } else if (wallstyle == ZCYLINDER) {
206 delxy = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1]);
207 if (delxy > 0.) {
208 delr = cylradius - delxy;
209
210 dx = -delr/delxy * x[i][0];
211 dy = -delr/delxy * x[i][1];
212
213 }
214 }
215
216 rsq = dx*dx + dy*dy + dz*dz;
217 if (rsq == 0.) continue; // center of the cylinder ... no repulsive force!
218
219 r = sqrt(rsq);
220 rinv = 1./r;
221
222 // repulsive penetration force
223 if (r <= r0) {
224
225 frac = r0*rinv;
226 frac2 = frac*frac;
227
228 fwall = D * frac2 *(frac2 - 1) * rinv;
229
230 f[i][0] += fwall * dx;
231 f[i][1] += fwall * dy;
232 f[i][2] += fwall * dz;
233
234 }
235
236 }
237 } // end loop nlocal
238 }
239
240 /* ---------------------------------------------------------------------- */
241
post_force_respa(int vflag,int ilevel,int iloop)242 void FixWallSph::post_force_respa(int vflag, int ilevel, int iloop)
243 {
244 if (ilevel == nlevels_respa-1) post_force(vflag);
245 }
246