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 Stefan Radl, TU Graz
36 Copyright 2016 - TU Graz
37 ------------------------------------------------------------------------- */
38
39 #include <cmath>
40 #include <stdlib.h>
41 #include <string.h>
42 #include "atom.h"
43 #include "atom_vec.h"
44 #include "comm.h"
45 #include "modify.h"
46 #include "force.h"
47 #include "memory.h"
48 #include "update.h"
49 #include "error.h"
50 #include "fix_mesh_surface.h"
51 #include "fix_neighlist_mesh.h"
52 #include "fix_multisphere.h"
53 #include "fix_property_atom.h"
54 #include "mpi_liggghts.h"
55 #include "math_extra_liggghts.h"
56 #include "fix_massflow_mesh_sieve.h"
57 #include "math_const.h"
58 #include "irregular.h"
59 #include "atom_vec_ellipsoid.h"
60 #include "random_park.h"
61
62 #include <algorithm>
63
64 using namespace LAMMPS_NS;
65 using namespace MathExtraLiggghts;
66 using namespace MathConst;
67 using namespace FixConst;
68
69 /* ---------------------------------------------------------------------- */
70
FixMassflowMeshSieve(LAMMPS * lmp,int narg,char ** arg)71 FixMassflowMeshSieve::FixMassflowMeshSieve(LAMMPS *lmp, int narg, char **arg) :
72 FixMassflowMesh(lmp, narg, arg),
73 sieveMultiSphereCanPass_(false),
74 verbose_(false), //developer to hard code here
75 sieveSize_(-1),
76 sieveSpacing_(-1),
77 sieveStiffness_(-1),
78 sieveDamping_(0),
79 fix_applyForce_(0),
80 random_(0)
81 {
82
83 // random number generator, seed is hardcoded
84 random_ = new RanPark(lmp,"15485863");
85
86 bool hasargs = true;
87
88 while(iarg_ < narg && hasargs)
89 {
90 hasargs = false;
91
92 if(strcmp(arg[iarg_],"sieveMultiSphereCanPass") == 0) {
93 sieveMultiSphereCanPass_ = true;
94 iarg_++;
95 hasargs = true;
96 }
97 else if ( strcmp(arg[iarg_],"sieveSize") == 0) {
98 if(narg < iarg_+2)
99 error->fix_error(FLERR,this,"not enough arguments for 'sieveSize'");
100 iarg_++;
101 sieveSize_ = atof(arg[iarg_++]);
102 hasargs = true;
103 }
104 else if ( strcmp(arg[iarg_],"sieveSpacing") == 0) {
105 if(narg < iarg_+2)
106 error->fix_error(FLERR,this,"not enough arguments for 'sieveSpacing'");
107 iarg_++;
108 sieveSpacing_ = atof(arg[iarg_++]);
109 hasargs = true;
110 }
111 else if ( strcmp(arg[iarg_],"sieveStiffness") == 0) {
112 if(narg < iarg_+2)
113 error->fix_error(FLERR,this,"not enough arguments for 'sieveStiffness'");
114 iarg_++;
115 sieveStiffness_ = atof(arg[iarg_++]);
116 hasargs = true;
117 }
118 else if ( strcmp(arg[iarg_],"sieveDamping") == 0) {
119 if(narg < iarg_+2)
120 error->fix_error(FLERR,this,"not enough arguments for 'sieveDamping'");
121 iarg_++;
122 sieveDamping_ = atof(arg[iarg_++]);
123 hasargs = true;
124 }
125 else
126 error->fix_error(FLERR,this,"unknown argument.");
127 }
128
129 if(sieveSize_<0 || sieveSize_> 1e5)
130 error->fix_error(FLERR,this,"'sieveSize' is not set, larger than 1e5, or less than zero.");
131
132 if(sieveSpacing_<0 || sieveSpacing_> 1e5)
133 error->fix_error(FLERR,this,"'sieveSpacing' is not set, larger than 1e5, or less than zero.");
134
135 if(sieveStiffness_<0 || sieveStiffness_> 1e10)
136 error->fix_error(FLERR,this,"'sieveStiffness' is not set, larger than 1e10, or less than zero.");
137 }
138
139 /* ---------------------------------------------------------------------- */
140
~FixMassflowMeshSieve()141 FixMassflowMeshSieve::~FixMassflowMeshSieve()
142 {
143 if (random_) delete random_;
144 }
145
146 /* ---------------------------------------------------------------------- */
post_create()147 void FixMassflowMeshSieve::post_create()
148 {
149
150 FixMassflowMesh::post_create();
151
152 // add per-particle counter for displacements
153 const char * fixarg[9];
154 sprintf(fixidApplyForce_,"massflowSieve_%s",id);
155 fixarg[0]=fixidApplyForce_;
156 fixarg[1]="all";
157 fixarg[2]="property/atom";
158 fixarg[3]=fixidApplyForce_;
159 fixarg[4]="scalar";
160 fixarg[5]="yes";
161 fixarg[6]="no";
162 fixarg[7]="no";
163 fixarg[8]="-1";
164 modify->add_fix(9,const_cast<char**>(fixarg));
165
166 fix_applyForce_ = static_cast<FixPropertyAtom*>(modify->find_fix_property(fixidApplyForce_,"property/atom","scalar",0,0,style));
167
168 }
169
170 /* ----------------------------------------------------------------------
171 initialize this fix
172 ------------------------------------------------------------------------- */
init()173 void FixMassflowMeshSieve::init()
174 {
175 FixMassflowMesh::init();
176 }
177
178 /* ---------------------------------------------------------------------- */
setmask()179 int FixMassflowMeshSieve::setmask()
180 {
181 int mask = FixMassflowMesh::setmask();
182 mask |= POST_FORCE;
183 mask |= PRE_EXCHANGE; //ensure pre_exchange always happens
184 return mask;
185 }
186
187 /* ----------------------------------------------------------------------
188 perform additinal force computations
189 ------------------------------------------------------------------------- */
post_force(int vflag)190 void FixMassflowMeshSieve::post_force(int vflag)
191 {
192 int nlocal = atom->nlocal;
193 double **x = atom->x;
194 double **v = atom->v;
195 double **f = atom->f;
196
197 double *radius = atom->radius;
198 int *mask = atom->mask;
199 double *counter = fix_counter_->vector_atom;
200 double *applyForce = fix_applyForce_->vector_atom;
201
202 double delta[3]={}, bary[3];
203 double deltan;
204
205 TriMesh *mesh = fix_mesh_->triMesh();
206 int nTriAll = mesh->sizeLocal() + mesh->sizeGhost();
207
208 // loop owned and ghost triangles
209 // count only if owned particle
210 for(int iTri = 0; iTri < nTriAll; iTri++)
211 {
212 const std::vector<int> & neighborList = fix_neighlist_->get_contact_list(iTri);
213 const int numneigh = neighborList.size();
214 double iTriDouble = double(iTri);
215
216 for(int iNeigh = 0; iNeigh < numneigh; iNeigh++)
217 {
218 const int iPart = neighborList[iNeigh];
219
220 // skip ghost particles
221 if(iPart >= nlocal)
222 continue;
223
224 // skip particles not in fix group
225 if (!(mask[iPart] & groupbit))
226 continue;
227
228 // in case of once_ == true, ignore everything which has been already counted
229 if(compDouble(counter[iPart],2.))
230 continue;
231
232 int barySign;
233 deltan = fix_mesh_->triMesh()->resolveTriSphereContactBary(iPart,iTri,radius[iPart],x[iPart],delta,bary,barySign);
234
235 if(deltan <= 0)
236 {
237 //Check if first contact, and throw coin to determine passage
238 if(applyForce[iPart]<0.0) //have first contact
239 {
240 double randNumber = random_->uniform();
241 if(!sieveMultiSphereCanPass_ && (fix_volumeweight_ms_ && fix_volumeweight_ms_->vector_atom[iPart] > 0)) //have an unwanted multisphere collision --> apply force in any case!
242 applyForce[iPart] = 1+iTriDouble;
243 else if( randNumber < sievePassProbability(radius[iPart]) )
244 {
245 applyForce[iPart] = 0;
246 continue;
247 }
248 else
249 applyForce[iPart] = 1+iTriDouble;
250
251 // printf("FixMassflowMeshSieve:: particle %d has an overlap with deltan: %g, delta: %g %g %g, radius: %g, applyForce: %g, randNumber: %g \n",
252 // iPart, deltan, delta[0],delta[1],delta[2],radius[iPart], applyForce[iPart],randNumber);
253
254 }
255 else if(compDouble(applyForce[iPart],0)) continue; //suppress force calculation
256
257 double contactNormal[3];
258 double invRadius = 1/radius[iPart];
259 contactNormal[0] = delta[0] * invRadius;
260 contactNormal[1] = delta[1] * invRadius;
261 contactNormal[2] = delta[2] * invRadius;
262
263 double
264 normalVel = v[iPart][0] * contactNormal[0]
265 + v[iPart][1] * contactNormal[1]
266 + v[iPart][2] * contactNormal[2];
267
268 const double Fn_damping = sieveDamping_ * normalVel;
269 const double Fn_contact =-sieveStiffness_* deltan; //deltan is NEGATIVE!
270 double Fn = Fn_damping + Fn_contact;
271 if(Fn<0.0) //avoid attractive forces during contact
272 Fn = 0.0;
273
274 f[iPart][0] -= Fn*contactNormal[0];
275 f[iPart][1] -= Fn*contactNormal[1];
276 f[iPart][2] -= Fn*contactNormal[2];
277 }
278 else //particle has no contact with this triangle now - but check if it had previously
279 {
280 if( compDouble(applyForce[iPart],1+iTriDouble) ) //previously a force was applied with this triangle, now release it!
281 applyForce[iPart] = -1;
282 //else --> have force with other triangle, or no interaction. Thus, no action for this triangle
283 }
284 }
285 }
286
287 }
288 /* ----------------------------------------------------------------------
289 perform some operations here if needed
290 done before exchange, borders, reneighbor
291 so that ghost atoms and neighbor lists will be correct
292 ------------------------------------------------------------------------- */
pre_exchange()293 void FixMassflowMeshSieve::pre_exchange()
294 {
295 FixMassflowMesh::pre_exchange();
296 }
297
sievePassProbability(double radius)298 double FixMassflowMeshSieve::sievePassProbability(double radius)
299 {
300 double dPass = sieveSize_ - 2.0 * radius;
301 if( dPass > 0.0 )
302 return dPass * dPass * PIOVERFOUR / (sieveSpacing_ * sieveSpacing_);
303 else
304 return 0;
305 };
306
307