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 no contributing author is listed, this file has been contributed
36 by the core developer)
37
38 Copyright 2015- DCS Computing GmbH, Linz
39 ------------------------------------------------------------------------- */
40
41 #include <cmath>
42 #include <string.h>
43 #include <stdlib.h>
44 #include "fix_base_liggghts.h"
45 #include "atom.h"
46 #include "error.h"
47 #include "force.h"
48 #include "group.h"
49 #include "update.h"
50 #include "respa.h"
51 #include "region.h"
52 #include "domain.h"
53 #include "modify.h"
54 #include "atom_vec.h"
55 #include "comm.h"
56 #include "fix_multisphere.h"
57 #include "mpi_liggghts.h"
58
59 using namespace LAMMPS_NS;
60 using namespace FixConst;
61
62 /* ---------------------------------------------------------------------- */
63
FixBaseLiggghts(LAMMPS * lmp,int narg,char ** arg)64 FixBaseLiggghts::FixBaseLiggghts(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, narg, arg),
65 support_respa_(false),
66 nlevels_respa_(0),
67 do_need_radius_(true),
68 do_need_mass_(true),
69 support_ms_(false),
70 fix_ms_(0),
71 ms_(0),
72 idregion_(0),
73 region_(0),
74 iregion_(-1)
75 {
76 }
77
78 /* ---------------------------------------------------------------------- */
79
~FixBaseLiggghts()80 FixBaseLiggghts::~FixBaseLiggghts()
81 {
82 if(idregion_) delete []idregion_;
83 }
84
85 /* ---------------------------------------------------------------------- */
86
process_region(char * regid)87 void FixBaseLiggghts::process_region(char *regid)
88 {
89 iregion_ = domain->find_region(regid);
90 if (iregion_ == -1)
91 error->fix_error(FLERR,this,"Region ID does not exist");
92 int n = strlen(regid) + 1;
93 idregion_ = new char[n];
94 strcpy(idregion_,regid);
95 region_ = domain->regions[iregion_];
96 }
97
98 /* ---------------------------------------------------------------------- */
99
init()100 void FixBaseLiggghts::init()
101 {
102 // error checks
103 if (do_need_radius_ && !atom->radius_flag)
104 error->fix_error(FLERR,this,"requires atom attribute radius (per-particle)");
105
106 // error checks
107 if (do_need_mass_ && !atom->rmass_flag)
108 error->fix_error(FLERR,this,"requires atom attribute mass (per-particle)");
109
110 // check validity of region
111 iregion_ = -1;
112 if (idregion_)
113 {
114 iregion_ = domain->find_region(idregion_);
115 if (iregion_ == -1)
116 error->fix_error(FLERR,this,"Region ID does not exist");
117 region_ = domain->regions[iregion_];
118 }
119
120 // multisphere support
121
122 fix_ms_ = static_cast<FixMultisphere*>(modify->find_fix_style("multisphere",0));
123 if(modify->n_fixes_style("multisphere") > 1)
124 error->fix_error(FLERR,this,"does not support more than one fix multisphere.");
125 if(fix_ms_)
126 ms_ = &(fix_ms_->data());
127 else
128 ms_ = 0;
129
130 if(!support_ms_ && ms_)
131 error->fix_error(FLERR,this,"does not support multi-sphere");
132
133 // currently no groups for MS
134 int igrp_all = group->find("all");
135 if(fix_ms_ && (igrp_all!= igroup))
136 error->fix_error(FLERR,this,"does only support fix group 'all' when multi-sphere particles present");
137
138 if (strstr(update->integrate_style,"respa"))
139 nlevels_respa_ = ((Respa *) update->integrate)->nlevels;
140 }
141
142 /* ---------------------------------------------------------------------- */
143
setup(int vflag)144 void FixBaseLiggghts::setup(int vflag)
145 {
146 if (strstr(update->integrate_style,"verlet"))
147 post_force(vflag);
148 else
149 {
150 if(!support_respa_)
151 error->fix_error(FLERR,this,"does nor support run_style respa");
152 ((Respa *) update->integrate)->copy_flevel_f(nlevels_respa_-1);
153 post_force_respa(vflag,nlevels_respa_-1,0);
154 ((Respa *) update->integrate)->copy_f_flevel(nlevels_respa_-1);
155 }
156 }
157
158 /* ---------------------------------------------------------------------- */
159
count_eligible(double & mass_counted,double & volume_counted,int & nparticles_counted)160 void FixBaseLiggghts::count_eligible(double &mass_counted,double &volume_counted, int &nparticles_counted)
161 {
162 int nlocal = atom->nlocal;
163 int nbody_local = ms_? ms_->n_body() : 0;
164 int *mask = atom->mask;
165 double **x = atom->x;
166 double *rmass = atom->rmass;
167 double *radius = atom->radius;
168 mass_counted = 0.;
169 volume_counted = 0.;
170 nparticles_counted = 0.;
171
172 double _4_pi_over_3 = 4.*M_PI/3.;
173
174 // spheres contribution
175
176 for(int ilocal = 0; ilocal < nlocal; ilocal++)
177 {
178 if (mask[ilocal] & groupbit)
179 {
180 // do not handle multi-sphere case
181 // skip if not in region
182 if (
183 (!fix_ms_ || fix_ms_->belongs_to(ilocal) < 0) &&
184 (!region_ || region_->match(x[ilocal][0],x[ilocal][1],x[ilocal][2]))
185 )
186 {
187 mass_counted += rmass[ilocal];
188 volume_counted += _4_pi_over_3*radius[ilocal]*radius[ilocal]*radius[ilocal];
189 nparticles_counted += 1;
190 }
191 }
192 }
193
194 // multi-spheres contribution
195
196 for(int ibody_local = 0; ibody_local < nbody_local ; ibody_local++)
197 {
198 double xcm[3];
199 ms_->xcm(xcm,ibody_local);
200
201 // skip if not in region; group not accounted for here
202 if (!region_ || region_->match(xcm[0],xcm[1],xcm[2]))
203 {
204 mass_counted += ms_->mass(ibody_local);
205 volume_counted += ms_->volume(ibody_local);
206 nparticles_counted += 1;
207 }
208 }
209
210 double vec_mpi[3];
211 vec_mpi[0] = mass_counted;
212 vec_mpi[1] = volume_counted;
213 vec_mpi[2] = static_cast<double>(nparticles_counted);
214 MPI_Sum_Vector(vec_mpi,3,world);
215 mass_counted = vec_mpi[0];
216 volume_counted = vec_mpi[1];
217 nparticles_counted = static_cast<int>(vec_mpi[2]);
218 }
219