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 <cmath>
43 #include <stdlib.h>
44 #include <string.h>
45 #include "atom.h"
46 #include "atom_vec.h"
47 #include "force.h"
48 #include "update.h"
49 #include "comm.h"
50 #include "modify.h"
51 #include "memory.h"
52 #include "error.h"
53 #include "group.h"
54 #include "neighbor.h"
55 #include "fix_property_global.h"
56
57 using namespace LAMMPS_NS;
58 using namespace FixConst;
59
60 #define EPSILON 0.001
61
62 /* ---------------------------------------------------------------------- */
63
FixPropertyGlobal(LAMMPS * lmp,int narg,char ** arg)64 FixPropertyGlobal::FixPropertyGlobal(LAMMPS *lmp, int narg, char **arg) :
65 Fix(lmp, narg, arg)
66 {
67 //Check args
68 if (narg < 6) error->all(FLERR,"Illegal fix property/global command, not enough arguments");
69
70 if(0 == strcmp(style,"custom_property/global"))
71 {
72 int len = strlen("property/global") + 1;
73 delete []style;
74 style = new char[len];
75 strcpy(style,"property/global");
76 }
77
78 //Read args
79 int n = strlen(arg[3]) + 1;
80 variablename = new char[n];
81 strcpy(variablename,arg[3]);
82 is_symmetric = false;
83 is_atomtype_bound = false;
84
85 if (strcmp(arg[4],"scalar") == 0)
86 data_style = FIXPROPERTY_GLOBAL_SCALAR;
87 else if (strcmp(arg[4],"vector") == 0)
88 data_style = FIXPROPERTY_GLOBAL_VECTOR;
89 else if (strcmp(arg[4],"peratomtype") == 0 || strcmp(arg[4],"atomtype") == 0)
90 {
91 data_style = FIXPROPERTY_GLOBAL_VECTOR;
92 is_atomtype_bound = true;
93 }
94 else if (strcmp(arg[4],"matrix") == 0)
95 data_style = FIXPROPERTY_GLOBAL_MATRIX;
96 else if (strcmp(arg[4],"peratomtypepair") == 0 || strcmp(arg[4],"atomtypepair") == 0)
97 {
98 data_style = FIXPROPERTY_GLOBAL_MATRIX;
99 is_symmetric = true;
100 is_atomtype_bound = true;
101 }
102 else error->fix_error(FLERR,this,"Unknown style. Valid styles are scalar, vector, atomtype/peratomtype, matrix, or atomtypepair/peratomtypepair");
103
104 int darg = 0;
105 if (data_style == FIXPROPERTY_GLOBAL_MATRIX) darg = 1;
106
107 //assign values
108 nvalues = narg - 5 - darg;
109 nvalues_new_array = 0;
110
111 values = (double*) memory->smalloc(nvalues*sizeof(double),"values");
112 values_recomputed = (double*) memory->smalloc(nvalues*sizeof(double),"values");
113
114 if(narg < 5+darg+nvalues) error->fix_error(FLERR,this,"not enough arguments");
115
116 for (int j = 0; j < nvalues; j++)
117 values[j] = force->numeric(FLERR,arg[5+darg+j]);
118
119 if (data_style == FIXPROPERTY_GLOBAL_SCALAR)
120 scalar_flag = 1;
121 else if (data_style==FIXPROPERTY_GLOBAL_VECTOR) {
122 vector_flag = 1;
123 size_vector = nvalues;
124 }
125 else if (data_style == FIXPROPERTY_GLOBAL_MATRIX) {
126 array_flag = 1;
127 size_array_cols = force->inumeric(FLERR,arg[5]);
128 if (fmod(static_cast<double>(nvalues),size_array_cols) != 0.)
129 error->fix_error(FLERR,this,"the number of default values must be a multiple of nCols.");
130 size_array_rows = static_cast<int>(static_cast<double>(nvalues)/size_array_cols);
131 }
132
133 extvector=0;
134
135 filename = 0;
136 grpname = 0;
137
138 //check if there is already a fix that tries to register a property with the same name
139
140 for (int ifix = 0; ifix < modify->nfix; ifix++)
141 if ((modify->fix[ifix]) && (strcmp(modify->fix[ifix]->style,style) == 0) && (strcmp(((FixPropertyGlobal*)(modify->fix[ifix]))->variablename,variablename)==0) )
142 error->fix_error(FLERR,this,"There is already a fix that registers a variable of the same name");
143
144 array = NULL;
145 array_recomputed = NULL;
146 if(data_style == FIXPROPERTY_GLOBAL_MATRIX)
147 {
148 array = (double**)memory->smalloc(size_array_rows*sizeof(double**),"FixPropGlob:array");
149 array_recomputed = (double**)memory->smalloc(size_array_rows*sizeof(double**),"FixPropGlob:array_recomputed");
150 for(int i = 0; i < size_array_rows; i++) array[i] = &values[i*size_array_cols];
151 for(int i = 0; i < size_array_rows; i++) array_recomputed[i] = &values_recomputed[i*size_array_cols];
152 }
153
154 // error check if matrix is symmetric (if required)
155 if(is_symmetric)
156 {
157 if(size_array_rows != size_array_cols)
158 error->fix_error(FLERR,this,"per-atomtype property matrix must be symmetric, i.e. N atom types "
159 "require you to define N columns and N rows with N*N total values");
160
161 int sflag = true;
162 for(int i = 0; i < size_array_rows; i++)
163 for(int j = 0; j < size_array_cols; j++)
164 if(array[i][j] != array[j][i])
165 sflag = false;
166
167 if(!lmp->wb)
168 {
169 if(!sflag)
170 error->fix_error(FLERR,this,"per-atomtype property matrix must be symmetric");
171 }
172 else
173 {
174 char errstr[512];
175 sprintf(errstr,"Property %s is required to be symmetric",variablename);
176 if(!sflag)
177 error->all(FLERR,errstr);
178 }
179 }
180 }
181
182 /* ---------------------------------------------------------------------- */
183
~FixPropertyGlobal()184 FixPropertyGlobal::~FixPropertyGlobal()
185 {
186 // delete locally stored arrays
187 delete[] variablename;
188
189 if(filename) delete[] filename;
190 if(grpname) delete[] grpname;
191
192 memory->sfree(values);
193 memory->sfree(values_recomputed);
194
195 if(array) memory->sfree(array);
196 if(array_recomputed) memory->sfree(array_recomputed);
197 }
198
199 /* ---------------------------------------------------------------------- */
200
pre_delete(bool unfixflag)201 void FixPropertyGlobal::pre_delete(bool unfixflag)
202 {
203 if(filename) write();
204 }
205
206 /* ---------------------------------------------------------------------- */
207
check_fix(const char * varname,const char * svmstyle,int len1,int len2,const char * caller,bool errflag)208 Fix* FixPropertyGlobal::check_fix(const char *varname,const char *svmstyle,int len1,int len2,const char *caller,bool errflag)
209 {
210 char errmsg[400];
211
212 if(strcmp(varname,variablename) == 0)
213 {
214 if(strcmp(svmstyle,"scalar") == 0) len1 = 1;
215
216 // check variable style
217 if(
218 (strcmp(svmstyle,"scalar") == 0 && data_style != FIXPROPERTY_GLOBAL_SCALAR) ||
219 ((strcmp(svmstyle,"vector") == 0 || strcmp(svmstyle,"peratomtype") == 0) && data_style != FIXPROPERTY_GLOBAL_VECTOR) ||
220 ((strcmp(svmstyle,"matrix") == 0 || strcmp(svmstyle,"peratomtypepair") == 0) && data_style != FIXPROPERTY_GLOBAL_MATRIX)
221 )
222 {
223 if(errflag)
224 {
225 sprintf(errmsg,"%s style required for fix property/global variable %s for usage with %s",svmstyle,varname,caller);
226 error->fix_error(FLERR,this,errmsg);
227 }
228 else return NULL;
229 }
230
231 // check length
232 if((nvalues < len1) && ((data_style != FIXPROPERTY_GLOBAL_MATRIX) || ((data_style == FIXPROPERTY_GLOBAL_MATRIX) && (size_array_cols < len2))))
233 {
234 if(errflag)
235 {
236 sprintf(errmsg,"Length not sufficient for variable %s for usage with %s",varname,caller);
237 error->fix_error(FLERR,this,errmsg);
238 }
239 else return NULL;
240 }
241
242 // success
243 return static_cast<Fix*>(this);
244 }
245 return NULL;
246 }
247
248 /* ---------------------------------------------------------------------- */
249
init()250 void FixPropertyGlobal::init()
251 {
252 me = comm->me;
253
254 char errmsg[300];
255 int ntypes = atom->ntypes;
256
257 if(FIXPROPERTY_GLOBAL_VECTOR == data_style && is_atomtype_bound && nvalues != ntypes)
258 {
259
260 sprintf(errmsg,"Fix property/global: Length not correct for variable %s, length should be equal to %d (= number of atom types)",
261 variablename,ntypes);
262 error->fix_error(FLERR,this,errmsg);
263 }
264 if(FIXPROPERTY_GLOBAL_MATRIX == data_style && is_atomtype_bound && nvalues != ntypes*ntypes)
265 {
266 sprintf(errmsg,"Fix property/global: Length not correct for variable %s, length should be equal to %d ( = number of atom types * number of atom types)",
267 variablename,ntypes*ntypes);
268 error->fix_error(FLERR,this,errmsg);
269 }
270 }
271
272 /* ---------------------------------------------------------------------- */
273
grow(int len1,int len2)274 void FixPropertyGlobal::grow(int len1, int len2)
275 {
276 if(data_style == FIXPROPERTY_GLOBAL_SCALAR) error->fix_error(FLERR,this,"Can not grow global property of type scalar");
277 else if(data_style == FIXPROPERTY_GLOBAL_VECTOR && len1 > nvalues)
278 {
279 memory->grow(values,len1,"FixPropertyGlobal:values");
280 }
281 else if(data_style == FIXPROPERTY_GLOBAL_MATRIX && len1*len2 > nvalues)
282 {
283 values = (double*) memory->srealloc(values,len1*len2*sizeof(double),"FixPropertyGlobal:values");
284 size_array_rows = len1;
285 size_array_cols = len2;
286 nvalues = len1*len2;
287 array = (double**)memory->srealloc(array,size_array_rows*sizeof(double**),"FixPropGlob:array");
288 for(int i = 0; i < size_array_rows; i++) array[i] = &values[i*size_array_cols];
289 }
290 }
291
292 /* ---------------------------------------------------------------------- */
293
compute_scalar()294 double FixPropertyGlobal::compute_scalar()
295 {
296 return values[0];
297 }
298
299 /* ---------------------------------------------------------------------- */
300
compute_vector(int i)301 double FixPropertyGlobal::compute_vector(int i)
302 {
303 if (i>(nvalues-1))error->fix_error(FLERR,this,"Trying to access vector, but index out of bounds");
304 return values[i];
305 }
306
vector_modify(int i,double val)307 void FixPropertyGlobal::vector_modify(int i,double val)
308 {
309 if (i>(nvalues-1))error->fix_error(FLERR,this,"Trying to access vector, but index out of bounds");
310 values_recomputed[i] = val;
311 }
312
compute_vector_modified(int i)313 double FixPropertyGlobal::compute_vector_modified(int i)
314 {
315 if (i>(nvalues-1))error->fix_error(FLERR,this,"Trying to access vector, but index out of bounds");
316 return values_recomputed[i];
317 }
318
319 /* ---------------------------------------------------------------------- */
320
compute_array(int i,int j)321 double FixPropertyGlobal::compute_array(int i, int j) //i is row, j is column
322 {
323 if (i>(size_array_rows-1))error->fix_error(FLERR,this,"Trying to access matrix, but row index out of bounds");
324 if (j>(size_array_cols-1))error->fix_error(FLERR,this,"Trying to access matrix, but column index out of bounds");
325
326 return array[i][j];
327 }
328
array_modify(int i,int j,double val)329 void FixPropertyGlobal::array_modify(int i, int j,double val) //i is row, j is column
330 {
331 if (i>(size_array_rows-1))error->fix_error(FLERR,this,"Trying to access matrix, but row index out of bounds");
332 if (j>(size_array_cols-1))error->fix_error(FLERR,this,"Trying to access matrix, but column index out of bounds");
333
334 array_recomputed[i][j] = val;
335 }
336
compute_array_modified(int i,int j)337 double FixPropertyGlobal::compute_array_modified(int i, int j) //i is row, j is column
338 {
339 if (i>(size_array_rows-1))error->fix_error(FLERR,this,"Trying to access matrix, but row index out of bounds");
340 if (j>(size_array_cols-1))error->fix_error(FLERR,this,"Trying to access matrix, but column index out of bounds");
341
342 return array_recomputed[i][j];
343 }
344
345 /* ---------------------------------------------------------------------- */
346
setmask()347 int FixPropertyGlobal::setmask()
348 {
349 int mask = 0;
350 return mask;
351 }
352
353 /* ----------------------------------------------------------------------
354 memory usage of local atom-based arrays
355 ------------------------------------------------------------------------- */
356
memory_usage()357 double FixPropertyGlobal::memory_usage()
358 {
359 double bytes = nvalues * sizeof(double);
360 return bytes;
361 }
362
363 /* ----------------------------------------------------------------------
364 memory usage of local atom-based arrays
365 ------------------------------------------------------------------------- */
366
new_array(int l1,int l2)367 void FixPropertyGlobal::new_array(int l1,int l2)
368 {
369
370 if (data_style == FIXPROPERTY_GLOBAL_MATRIX) error->fix_error(FLERR,this,"Can not allocate extra array for matrix style");
371 array_flag = 1;
372 size_array_rows = l1;
373 size_array_cols = l2;
374 nvalues_new_array = l1*l2;
375
376 memory->create(array,size_array_rows,size_array_cols,"FixPropGlob:array");
377 memory->create(array_recomputed,size_array_rows,size_array_cols,"FixPropGlob:array_recomputed");
378 }
379
380 /* ----------------------------------------------------------------------
381 write out command
382 ------------------------------------------------------------------------- */
383
write()384 void FixPropertyGlobal::write()
385 {
386
387 if(0 != me)
388 return;
389
390 FILE *file = fopen(filename,"w");
391
392 if(!file)
393 error->one(FLERR,"Fix property/global cannot open file");
394
395 // fix id group style variablename
396 fprintf(file,"fix %s %s %s %s ",id,grpname,style,variablename);
397
398 // datatype
399 const char *datatyp;
400 if(0 == data_style) datatyp = "scalar";
401 if(1 == data_style) datatyp = "vector";
402 if(2 == data_style && is_symmetric) datatyp = "atomtypepair";
403 else if(2 == data_style) datatyp = "matrix";
404 fprintf(file,"%s ",datatyp);
405
406 // size_array_cols if required
407 if(2 == data_style) fprintf(file,"%d ",size_array_cols);
408
409 // values
410 for(int i = 0; i < nvalues; i++)
411 fprintf(file,"%f ",values[i]);
412
413 fprintf(file,"\n");
414 fclose(file);
415 }
416
417 /* ---------------------------------------------------------------------- */
418
modify_param(int narg,char ** arg)419 int FixPropertyGlobal::modify_param(int narg, char **arg)
420 {
421 if (strcmp(arg[0],"file") == 0) {
422 if (narg < 2) error->fix_error(FLERR,this,"not enough arguments for fix_modify 'file'");
423
424 filename = new char[strlen(arg[1])+1];
425 strcpy(filename,arg[1]);
426 grpname = new char[strlen(group->names[igroup])+1];
427 strcpy(grpname,group->names[igroup]);
428 return 2;
429 }
430
431 return 0;
432 }
433