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