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