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 (if not contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2012- DCS Computing GmbH, Linz
39 Copyright 2009-2012 JKU Linz
40 ------------------------------------------------------------------------- */
41
42 #include "particleToInsert.h"
43 #include <cmath>
44 #include "error.h"
45 #include "update.h"
46 #include "domain.h"
47 #include "atom.h"
48 #include "atom_vec.h"
49 #include "fix_property_atom.h"
50 #include "vector_liggghts.h"
51 #include "math_extra_liggghts.h"
52 #include "modify.h"
53
54 using namespace LAMMPS_NS;
55
ParticleToInsert(LAMMPS * lmp,int ns,FixPropertyAtom * const _fix_release)56 ParticleToInsert::ParticleToInsert(LAMMPS* lmp, int ns, FixPropertyAtom * const _fix_release) :
57 Pointers(lmp),
58 id_ins(-1),
59 fix_release(_fix_release),
60 fix_property(NULL),
61 n_fix_property(0),
62 fix_property_nentry(NULL),
63 fix_property_value(NULL),
64 fix_template_(NULL)
65 {
66 groupbit = 0;
67
68 distorder = -1;
69
70 nparticles = ns;
71
72 memory->create(x_ins,nparticles,3,"x_ins");
73 radius_ins = new double[nparticles];
74
75 atom_type_vector = new int[nparticles];
76 atom_type_vector_flag = false;
77
78 if (fix_release && fix_release->get_nvalues() <= 14)
79 error->all(FLERR, "Invalid fix_release, more than 14 entries required");
80 }
81
82 /* ---------------------------------------------------------------------- */
83
~ParticleToInsert()84 ParticleToInsert::~ParticleToInsert()
85 {
86 memory->destroy(x_ins);
87 delete []radius_ins;
88 delete []atom_type_vector;
89 if (fix_property_value)
90 {
91 for (int i = 0; i < n_fix_property; i++)
92 delete [] fix_property_value[i];
93 delete [] fix_property_value;
94 }
95 if (fix_property)
96 delete [] fix_property;
97 }
98
99 /* ---------------------------------------------------------------------- */
100
insert()101 int ParticleToInsert::insert()
102 {
103 // perform the actual insertion
104 // add particles, set coordinate and radius
105 // set group mask to "all" plus fix groups
106
107 int inserted = 0;
108 int nfix = modify->nfix;
109 Fix **fix = modify->fix;
110
111 for(int i = 0; i < nparticles; i++)
112 {
113
114 //if (domain->is_in_extended_subdomain(x_ins[i]))
115 //{
116
117 inserted++;
118 if(atom_type_vector_flag)
119 atom->avec->create_atom(atom_type_vector[i],x_ins[i]);
120 else
121 atom->avec->create_atom(atom_type,x_ins[i]);
122 int m = atom->nlocal - 1;
123 atom->mask[m] = 1 | groupbit;
124 vectorCopy3D(v_ins,atom->v[m]);
125 vectorCopy3D(omega_ins,atom->omega[m]);
126 atom->radius[m] = radius_ins[i];
127 atom->density[m] = density_ins;
128
129 atom->rmass[m] = (1==nparticles)? (mass_ins) : (4.18879020479/*4//3*pi*/*radius_ins[i]*radius_ins[i]*radius_ins[i]*density_ins);
130
131 //pre_set_arrays() called via FixParticleDistribution
132 for (int j = 0; j < nfix; j++)
133 if (fix[j]->create_attribute) fix[j]->set_arrays(m);
134
135 // apply fix property setting coming from fix insert
136 // this overrides the set_arrays call above
137 if(fix_property)
138 {
139 for (int j = 0; j < n_fix_property; j++)
140 {
141 if (fix_property_nentry[j] == 1)
142 fix_property[j]->vector_atom[m] = fix_property_value[j][0];
143 else
144 {
145 for (int k = 0; k < fix_property_nentry[j]; k++)
146 fix_property[j]->array_atom[m][k] = fix_property_value[j][k];
147 }
148 }
149 }
150 if (fix_template_)
151 fix_template_->vector_atom[m] = (double)distorder;
152 if (fix_release)
153 fix_release->array_atom[m][14] = (double) id_ins;
154 //}
155 }
156
157 return inserted;
158 }
159
160 /* ---------------------------------------------------------------------- */
161
162 /*
163 int ParticleToInsert::check_near_set_x_v_omega(double *x,double *v, double *omega, double *quat, double **xnear, int &nnear)
164 {
165 if(nparticles > 1)
166 return check_near_set_x_v_omega_ms(x,v, omega,quat,xnear,nnear);
167
168 // check sphere against all others in xnear
169 // if no overlap add to xnear
170 double del[3], rsq, radsum;
171
172 vectorCopy3D(x,x_ins[0]);
173
174 for(int i = 0; i < nnear; i++)
175 {
176 vectorSubtract3D(x_ins[0],xnear[i],del);
177 rsq = vectorMag3DSquared(del);
178
179 /*
180 radsum = radius_ins[0] + xnear[i][3];
181
182 // no success in overlap
183 if (rsq <= radsum*radsum) return 0;
184 }
185
186 // no overlap with any other - success
187
188 vectorCopy3D(v,v_ins);
189 vectorCopy3D(omega,omega_ins);
190
191 // add to xnear
192 vectorCopy3D(x_ins[0],xnear[nnear]);
193 xnear[nnear][3] = radius_ins[0];
194 nnear++;
195
196 return 1;
197 }*/
198
199 /* ---------------------------------------------------------------------- */
200
201 /*
202 int ParticleToInsert::check_near_set_x_v_omega_ms(double *x,double *v, double *omega, double *quat, double **xnear, int &nnear)
203 {
204 // x is position where insertion should take place
205 // v and omega are the velocity and omega for the newly inserted particles
206 double rel[3],xins_j_try[3];
207 double del[3], rsq, radsum;
208
209 // check insertion position, take quat into account
210 // relative position of spheres to each other already stored at this point
211 // check sphere against all others in xnear
212 for(int j = 0; j < nspheres; j++)
213 {
214 // take orientation into account; x_bound_ins is in the global coordinate system
215 // calculate xins_j_try for every sphere and check if would work
216 vectorSubtract3D(x_ins[j],x_bound_ins,rel);
217 MathExtraLiggghts::vec_quat_rotate(rel,quat);
218 vectorAdd3D(rel,x,xins_j_try);
219
220 for(int i = 0; i < nnear; i++)
221 {
222 vectorSubtract3D(xins_j_try,xnear[i],del);
223 rsq = vectorMag3DSquared(del);
224 radsum = radius_ins[j] + xnear[i][3];
225
226 // no success in overlap
227 if (rsq <= radsum*radsum)
228 return 0;
229 }
230 }
231
232 // no overlap with any other - success
233 // set x_ins, v_ins and omega_ins
234 for(int j = 0; j < nspheres; j++)
235 {
236 vectorSubtract3D(x_ins[j],x_bound_ins,rel);
237 MathExtraLiggghts::vec_quat_rotate(rel,quat);
238 vectorAdd3D(rel,x,x_ins[j]);
239 }
240 vectorCopy3D(v,v_ins);
241 vectorCopy3D(omega,omega_ins);
242
243 // add to xnear for future checks
244 for(int j = 0; j < nspheres; j++)
245 {
246 vectorCopy3D(x_ins[j],xnear[nnear]);
247 xnear[nnear][3] = radius_ins[j];
248 nnear++;
249 }
250
251 return nparticles;
252 }*/
253
254 /* ---------------------------------------------------------------------- */
255
check_near_set_x_v_omega(double * x,double * v,double * omega,double * quat,RegionNeighborList<interpolate_no> & neighList)256 int ParticleToInsert::check_near_set_x_v_omega(double *x,double *v, double *omega, double *quat, RegionNeighborList<interpolate_no> & neighList)
257 {
258 if(nparticles > 1)
259 return check_near_set_x_v_omega_ms(x,v, omega,quat,neighList);
260
261 vectorCopy3D(x,x_ins[0]);
262
263 if(neighList.hasOverlap(x_ins[0], radius_ins[0])) {
264 return 0;
265 }
266
267 // no overlap with any other - success
268
269 vectorCopy3D(v,v_ins);
270 vectorCopy3D(omega,omega_ins);
271
272 neighList.insert(x_ins[0], radius_ins[0]);
273
274 return 1;
275 }
276
277 /* ---------------------------------------------------------------------- */
278
check_near_set_x_v_omega_ms(double * x,double * v,double * omega,double * quat,RegionNeighborList<interpolate_no> & neighList)279 int ParticleToInsert::check_near_set_x_v_omega_ms(double *x,double *v, double *omega, double *quat, RegionNeighborList<interpolate_no> & neighList)
280 {
281 // x is position where insertion should take place
282 // v and omega are the velocity and omega for the newly inserted particles
283 double rel[3],xins_j_try[3];
284 //double del[3], rsq, radsum;
285
286 // check insertion position, take quat into account
287 // relative position of spheres to each other already stored at this point
288 // check sphere against all others in xnear
289 for(int j = 0; j < nparticles; j++)
290 {
291 // take orientation into account; x_bound_ins is in the global coordinate system
292 // calculate xins_j_try for every sphere and check if would work
293 vectorSubtract3D(x_ins[j],x_bound_ins,rel);
294 MathExtraLiggghts::vec_quat_rotate(rel,quat);
295 vectorAdd3D(rel,x,xins_j_try);
296
297 if(neighList.hasOverlap(xins_j_try, radius_ins[j])) {
298 return 0;
299 }
300 }
301
302 // no overlap with any other - success
303 // set x_ins, v_ins and omega_ins
304 for(int j = 0; j < nparticles; j++)
305 {
306 vectorSubtract3D(x_ins[j],x_bound_ins,rel);
307 MathExtraLiggghts::vec_quat_rotate(rel,quat);
308 vectorAdd3D(rel,x,x_ins[j]);
309 }
310 vectorCopy3D(v,v_ins);
311 vectorCopy3D(omega,omega_ins);
312
313 // add to xnear for future checks
314 for(int j = 0; j < nparticles; j++)
315 {
316 neighList.insert(x_ins[j], radius_ins[j]);
317 }
318
319 return nparticles;
320 }
321
322 /* ---------------------------------------------------------------------- */
323
set_x_v_omega(double * x,double * v,double * omega,double * quat)324 int ParticleToInsert::set_x_v_omega(double *x, double *v, double *omega, double *quat)
325 {
326 double rel[3];
327
328 // x is position where insertion should take place
329 // v and omega are the velocity and omega for the newly inserted particles
330
331 // add insertion position
332 // relative position of spheres to each other already stored at this point
333 // also take quat into account
334 for(int j = 0; j < nparticles; j++)
335 {
336 // if only one sphere, then x_bound = x_ins and there is
337 // no relevant orientation
338 if(1 == nparticles)
339 vectorAdd3D(x_ins[j],x,x_ins[j]);
340 // if > 1 sphere, take orientation into account
341 // x_bound_ins is in the global coordinate system
342 else
343 {
344 vectorSubtract3D(x_ins[j],x_bound_ins,rel);
345 MathExtraLiggghts::vec_quat_rotate(rel,quat);
346 vectorAdd3D(rel,x,x_ins[j]);
347 }
348 }
349
350 // set velocity and omega
351 vectorCopy3D(v,v_ins);
352 vectorCopy3D(omega,omega_ins);
353
354 return nparticles;
355 }
356
357 /* ---------------------------------------------------------------------- */
358
scale_pti(double r_scale)359 void ParticleToInsert::scale_pti(double r_scale)
360 {
361 double r_scale3 = r_scale*r_scale*r_scale;
362
363 for(int i = 0; i < nparticles; i++)
364 {
365 radius_ins[i] *= r_scale;
366 vectorScalarMult3D(x_ins[i],r_scale);
367 }
368
369 volume_ins *= r_scale3;
370 mass_ins *= r_scale3;
371
372 r_bound_ins *= r_scale;
373
374 vectorScalarMult3D(x_bound_ins,r_scale);
375 }
376