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 This file is from LAMMPS
36 LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
37 http://lammps.sandia.gov, Sandia National Laboratories
38 Steve Plimpton, sjplimp@sandia.gov
39
40 Copyright (2003) Sandia Corporation. Under the terms of Contract
41 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
42 certain rights in this software. This software is distributed under
43 the GNU General Public License.
44 ------------------------------------------------------------------------- */
45
46 #include <stdlib.h>
47 #include <string.h>
48 #include "compute_slice.h"
49 #include "update.h"
50 #include "modify.h"
51 #include "fix.h"
52 #include "group.h"
53 #include "memory.h"
54 #include "error.h"
55 #include "force.h"
56
57 using namespace LAMMPS_NS;
58
59 enum{COMPUTE,FIX};
60
61 /* ---------------------------------------------------------------------- */
62
ComputeSlice(LAMMPS * lmp,int & iarg,int narg,char ** arg)63 ComputeSlice::ComputeSlice(LAMMPS *lmp, int &iarg, int narg, char **arg) :
64 Compute(lmp, iarg, narg, arg)
65 {
66 if (narg < iarg+4) error->all(FLERR,"Illegal compute slice command");
67
68 MPI_Comm_rank(world,&me);
69
70 nstart = force->inumeric(FLERR,arg[iarg++]);
71 nstop = force->inumeric(FLERR,arg[iarg++]);
72 nskip = force->inumeric(FLERR,arg[iarg++]);
73
74 if (nstart < 1 || nstop < nstart || nskip < 1)
75 error->all(FLERR,"Illegal compute slice command");
76
77 // parse remaining values until one isn't recognized
78
79 which = new int[narg-iarg];
80 argindex = new int[narg-iarg];
81 ids = new char*[narg-iarg];
82 value2index = new int[narg-iarg];
83 nvalues = 0;
84
85 for (; iarg < narg; iarg++) {
86 if (strncmp(arg[iarg],"c_",2) == 0 ||
87 strncmp(arg[iarg],"f_",2) == 0) {
88 if (arg[iarg][0] == 'c') which[nvalues] = COMPUTE;
89 else if (arg[iarg][0] == 'f') which[nvalues] = FIX;
90
91 int n = strlen(arg[iarg]);
92 char *suffix = new char[n];
93 strcpy(suffix,&arg[iarg][2]);
94
95 char *ptr = strchr(suffix,'[');
96 if (ptr) {
97 if (suffix[strlen(suffix)-1] != ']')
98 error->all(FLERR,"Illegal compute slice command");
99 argindex[nvalues] = atoi(ptr+1);
100 *ptr = '\0';
101 } else argindex[nvalues] = 0;
102
103 n = strlen(suffix) + 1;
104 ids[nvalues] = new char[n];
105 strcpy(ids[nvalues],suffix);
106 nvalues++;
107 delete [] suffix;
108
109 } else error->all(FLERR,"Illegal compute slice command");
110 }
111
112 // setup and error check
113
114 for (int i = 0; i < nvalues; i++) {
115 if (which[i] == COMPUTE) {
116 int icompute = modify->find_compute(ids[i]);
117 if (icompute < 0)
118 error->all(FLERR,"Compute ID for compute slice does not exist");
119 if (modify->compute[icompute]->vector_flag) {
120 if (argindex[i])
121 error->all(FLERR,"Compute slice compute does not calculate a global array");
122 if (nstop > modify->compute[icompute]->size_vector)
123 error->all(FLERR,"Compute slice compute vector is accessed out-of-range");
124 } else if (modify->compute[icompute]->array_flag) {
125 if (argindex[i] == 0)
126 error->all(FLERR,"Compute slice compute does not calculate a global vector");
127 if (argindex[i] > modify->compute[icompute]->size_array_cols)
128 error->all(FLERR,"Compute slice compute array is accessed out-of-range");
129 if (nstop > modify->compute[icompute]->size_array_rows)
130 error->all(FLERR,"Compute slice compute array is accessed out-of-range");
131 } else error->all(FLERR,"Compute slice compute does not calculate "
132 "global vector or array");
133 } else if (which[i] == FIX) {
134 int ifix = modify->find_fix(ids[i]);
135 if (ifix < 0)
136 error->all(FLERR,"Fix ID for compute slice does not exist");
137 if (modify->fix[ifix]->vector_flag) {
138 if (argindex[i])
139 error->all(FLERR,"Compute slice fix does not calculate a global array");
140 if (nstop > modify->fix[ifix]->size_vector)
141 error->all(FLERR,"Compute slice fix vector is accessed out-of-range");
142 } else if (modify->fix[ifix]->array_flag) {
143 if (argindex[i] == 0)
144 error->all(FLERR,"Compute slice fix does not calculate a global vector");
145 if (argindex[i] > modify->fix[ifix]->size_array_cols)
146 error->all(FLERR,"Compute slice fix array is accessed out-of-range");
147 if (nstop > modify->fix[ifix]->size_array_rows)
148 error->all(FLERR,"Compute slice fix array is accessed out-of-range");
149 } else error->all(FLERR,"Compute slice fix does not calculate "
150 "global vector or array");
151 }
152 }
153
154 // this compute produces either a vector or array
155 // for vector, set intensive/extensive to mirror input values
156 // for array, set intensive if all input values are intensive, else extensive
157
158 vector = NULL;
159 array = NULL;
160 extlist = NULL;
161
162 if (nvalues == 1) {
163 vector_flag = 1;
164 size_vector = (nstop-nstart) / nskip;
165 memory->create(vector,size_vector,"slice:vector");
166
167 if (which[0] == COMPUTE) {
168 int icompute = modify->find_compute(ids[0]);
169 if (argindex[0] == 0) {
170 extvector = modify->compute[icompute]->extvector;
171 if (modify->compute[icompute]->extvector == -1) {
172 extlist = new int[size_vector];
173 int j = 0;
174 for (int i = nstart; i < nstop; i += nskip)
175 extlist[j++] = modify->compute[icompute]->extlist[i-1];
176 }
177 } else extvector = modify->compute[icompute]->extarray;
178 } else if (which[0] == FIX) {
179 int ifix = modify->find_fix(ids[0]);
180 if (argindex[0] == 0) {
181 extvector = modify->fix[ifix]->extvector;
182 if (modify->fix[ifix]->extvector == -1) {
183 extlist = new int[size_vector];
184 int j = 0;
185 for (int i = nstart; i < nstop; i += nskip)
186 extlist[j++] = modify->fix[ifix]->extlist[i-1];
187 }
188 } else extvector = modify->fix[ifix]->extarray;
189 }
190
191 } else {
192 array_flag = 1;
193 size_array_rows = (nstop-nstart) / nskip;
194 size_array_cols = nvalues;
195 memory->create(array,size_array_rows,size_array_cols,"slice:array");
196
197 extarray = 0;
198 for (int i = 0; i < nvalues; i++) {
199 if (which[i] == COMPUTE) {
200 int icompute = modify->find_compute(ids[i]);
201 if (argindex[i] == 0) {
202 if (modify->compute[icompute]->extvector == 1) extarray = 1;
203 if (modify->compute[icompute]->extvector == -1) {
204 for (int j = 0; j < modify->compute[icompute]->size_vector; j++)
205 if (modify->compute[icompute]->extlist[j]) extarray = 1;
206 }
207 } else {
208 if (modify->compute[icompute]->extarray) extarray = 1;
209 }
210 } else if (which[i] == FIX) {
211 int ifix = modify->find_fix(ids[i]);
212 if (argindex[i] == 0) {
213 if (modify->fix[ifix]->extvector == 1) extarray = 1;
214 if (modify->fix[ifix]->extvector == -1) {
215 for (int j = 0; j < modify->fix[ifix]->size_vector; j++)
216 if (modify->fix[ifix]->extlist[j]) extarray = 1;
217 }
218 } else {
219 if (modify->fix[ifix]->extarray) extarray = 1;
220 }
221 }
222 }
223 }
224 }
225
226 /* ---------------------------------------------------------------------- */
227
~ComputeSlice()228 ComputeSlice::~ComputeSlice()
229 {
230 delete [] which;
231 delete [] argindex;
232 for (int m = 0; m < nvalues; m++) delete [] ids[m];
233 delete [] ids;
234 delete [] value2index;
235
236 memory->destroy(vector);
237 memory->destroy(array);
238 }
239
240 /* ---------------------------------------------------------------------- */
241
init()242 void ComputeSlice::init()
243 {
244 // set indices and check validity of all computes,fixes
245
246 for (int m = 0; m < nvalues; m++) {
247 if (which[m] == COMPUTE) {
248 int icompute = modify->find_compute(ids[m]);
249 if (icompute < 0)
250 error->all(FLERR,"Compute ID for compute slice does not exist");
251 value2index[m] = icompute;
252 } else if (which[m] == FIX) {
253 int ifix = modify->find_fix(ids[m]);
254 if (ifix < 0)
255 error->all(FLERR,"Fix ID for compute slice does not exist");
256 value2index[m] = ifix;
257 }
258 }
259 }
260
261 /* ---------------------------------------------------------------------- */
262
compute_vector()263 void ComputeSlice::compute_vector()
264 {
265 invoked_vector = update->ntimestep;
266
267 extract_one(0,vector,1);
268 }
269
270 /* ---------------------------------------------------------------------- */
271
compute_array()272 void ComputeSlice::compute_array()
273 {
274 invoked_array = update->ntimestep;
275
276 for (int m = 0; m < nvalues; m++)
277 extract_one(0,&array[m][0],nvalues);
278 }
279
280 /* ----------------------------------------------------------------------
281 calculate sliced value for one input M and return it in vec
282 vec may be array so that returned values are with stride
283 ------------------------------------------------------------------------- */
284
extract_one(int m,double * vec,int stride)285 void ComputeSlice::extract_one(int m, double *vec, int stride)
286 {
287 int i,j;
288
289 // invoke the appropriate compute if needed
290
291 if (which[m] == COMPUTE) {
292 Compute *compute = modify->compute[value2index[m]];
293
294 if (argindex[m] == 0) {
295 if (!(compute->invoked_flag & INVOKED_VECTOR)) {
296 compute->compute_vector();
297 compute->invoked_flag |= INVOKED_VECTOR;
298 }
299 double *cvector = compute->vector;
300 j = 0;
301 for (i = nstart; i < nstop; i += nskip) {
302 vec[j] = cvector[i-1];
303 j += stride;
304 }
305
306 } else {
307 if (!(compute->invoked_flag & INVOKED_ARRAY)) {
308 compute->compute_array();
309 compute->invoked_flag |= INVOKED_ARRAY;
310 }
311 double **carray = compute->array;
312 int icol = argindex[m]-1;
313 j = 0;
314 for (i = nstart; i < nstop; i += nskip) {
315 vec[j] = carray[i-1][icol];
316 j += stride;
317 }
318 }
319
320 // access fix fields, check if fix frequency is a match
321
322 } else if (which[m] == FIX) {
323 if (update->ntimestep % modify->fix[value2index[m]]->global_freq)
324 error->all(FLERR,"Fix used in compute slice not computed at compatible time");
325 Fix *fix = modify->fix[value2index[m]];
326
327 if (argindex[m] == 0) {
328 j = 0;
329 for (i = nstart; i < nstop; i += nskip) {
330 vec[j] = fix->compute_vector(i-1);
331 j += stride;
332 }
333 } else {
334 int icol = argindex[m]-1;
335 j = 0;
336 for (i = nstart; i < nstop; i += nskip) {
337 vec[j] = fix->compute_array(i-1,icol);
338 j += stride;
339 }
340 }
341 }
342 }
343