NDDEM
Parameters.h
Go to the documentation of this file.
1 #ifndef PARAMETERS_H
2 #define PARAMETERS_H
3 
4 //using json = nlohmann::json;
5 using json = nddem::json ;
6 
9 //enum DataValue {radius, mass, Imom, pos, vel, omega, id1, id2, pospq, fpq, mpq, mqp} ;
10 
11 class Param {
12 public:
13 
14  int skipT=0 ;
15  int maxT =-1;
16  vector <string> flags = {} ;
17  int dim=-1 ;
18  vector <int> boxes ;
19  vector <vector <double> > boundaries ;
20  string save="" ;
21  vector<string> saveformat ;
23 
25 
27  double windowsize ;
28  vector<bool> periodicity ;
29  vector<double>delta ;
30 
31  std::map <std::string, bool> requiredfieldset = {{"dimension", false}, {"window size", false}, {"fields", false}, {"max time", false}, {"boundaries", false}, {"boxes", false}} ;
32 
33  //----- Files
34  struct File {
37  std::map <std::string, std::string> mapping ;
39  //~File() {if (reader != nullptr) {reader->close() ; delete (reader) ; }}
40  } ;
41  std::vector <File> files ;
42 
43  //----- Extra fields
44  struct ExtaField {
45  string name ;
47  FieldType type ;
48  std::optional<std::string> mapping ;
49  int datalocation ;
50  } ;
51  vector<ExtaField> extrafields ;
52 
53  int curts=0 ;
54  bool tsread = false ;
55 
56  void from_json (json & j) ;
57  int read_timestep (int ts, bool particleonly=false) ;
58  int set_data (Data & cgdata) ;
59  int get_num_particles () ;
60  int get_num_contacts () ;
61  double * get_data(DataValue datavalue, int dd=0, std::string name="") ;
62  void post_init () ;
63 
64 
65 private :
66  int identify_max_level() ;
67  Windows identify_window(std::string windowname) ;
68  void process_extrafields (json &j) ;
69  void process_file (json &j) ;
70 } ;
71 
72 //=================================================================================================
73 void Param::from_json(json &j)
74 {
75  for (auto v : j.items())
76  {
77  if (v.key() == "file") {process_file(v.value()) ; }
78  else if (v.key() == "savefile") save = v.value().get<string>() ;
79  else if (v.key() == "saveformat")
80  {
81  try {saveformat = v.value().get<vector<string>>();}
82  catch(...) {string a = v.value().get<string>(); saveformat.push_back(a) ; }
83  }
84  else if (v.key() == "window size")
85  { windowsize = v.value().get<double>(); requiredfieldset[v.key()] = true ; }
86  else if (v.key() == "skip") skipT = v.value().get<int>() ;
87  else if (v.key() == "max time") { maxT = v.value().get<int>() ; requiredfieldset[v.key()] = true ;}
88  else if (v.key() == "fields") { flags = v.value().get<vector<string>>(); requiredfieldset[v.key()] = true ; }
89  else if (v.key() == "boxes") {boxes = v.value().get<decltype(boxes)>(); requiredfieldset[v.key()] = true ; }
90  else if (v.key() == "boundaries") {boundaries = v.value().get<decltype(boundaries)>(); requiredfieldset[v.key()] = true ; }
91  else if (v.key() == "time average")
92  {
93  if (v.value().get<string>() == "None") timeaverage = AverageType::None ;
94  else if (v.value().get<string>() == "Final") timeaverage = AverageType::Final ;
95  else if (v.value().get<string>() == "Intermediate") timeaverage = AverageType::Intermediate ;
96  else if (v.value().get<string>() == "Intermediate and Final" || v.value().get<string>() =="Final and Intermediate") timeaverage = AverageType::Both ;
97  else
98  printf("WARN: unknown average type (allowed: 'None', 'Final', 'Intermediate', 'Intermediate and Final'). Case sensitive\n") ;
99  }
100  else if (v.key() == "window") {window = identify_window(v.value().get<string>()); }
101  else if (v.key() == "periodicity") {periodicity = v.value().get<decltype(periodicity)>();}
102  else if (v.key() == "density") { default_density=v.value().get<double>() ; }
103  else if (v.key() == "radius") { default_radius=v.value().get<double>() ; }
104  else if (v.key() == "diameter") { default_radius=v.value().get<double>()/2. ; }
105  else if (v.key() == "extra fields") { process_extrafields(v.value()) ; }
106  else if (v.key() == "dimension") {dim = v.value().get<int>() ; requiredfieldset[v.key()] = true ; }
107  else printf("Unknown json key: %s.\n", v.key().c_str()) ;
108  }
109 
110 
111 /*
112  post_init() ;
113 */
114 }
115 //-----------------------------------------------
116 int Param::read_timestep (int ts, bool particleonly)
117 {
118  if (curts == ts && tsread) return 0 ; //Already read and all
119  printf("%d ", ts) ;fflush(stdout) ;
120  for (auto & v: files)
121  {
122  v.reader->read_timestep(ts) ;
123  }
124  return 1 ;
125 }
126 //-------------------------------
127 void Param::post_init()
128 {
129  if (files.size() ==0) { printf("ERR: a file should be defined.\n"); return ; }
130  if (dim == -1 || !requiredfieldset["dimension"] ) dim=files[0].reader->get_dimension() ;
131  if (!requiredfieldset["window size"]) { printf("ERR: 'window size' is required and cannot be built. \n") ; return ; }
132  if (!requiredfieldset["fields"]) { printf("ERR: 'fields' is required and cannot be built. \n") ; return ; }
133  if (!requiredfieldset["boxes"]) { printf("ERR: the number of boxes in each dimension is required and cannot be built.\n") ; return ; }
134  if (!requiredfieldset["max time"])
135  {
136  maxT = files[0].reader->get_numts() ;
137  printf("//%d %d//", maxT, skipT) ;
138  if (maxT==-1)
139  {
140  printf("ERR: Cannot find the total number of timestep from the file.\n") ;
141  return ;
142  }
143  maxT -= skipT ;
144  printf("//%d//", maxT) ;
145  if (maxT<0) maxT=0 ;
146  }
147  if (!requiredfieldset["boundaries"]) boundaries= files[0].reader->get_bounds() ;
148 
149  if (boundaries.size() != 2 || boundaries[0].size() != static_cast<unsigned int>(dim) || boundaries[1].size()!=static_cast<unsigned int>(dim))
150  printf("ERR: dimension of the boundaries is not consistent (should be '2 x dimension')\n") ;
151  if (boxes.size() != static_cast<unsigned int>(dim))
152  printf("ERR: dimension of the boxes is not consistent (should be 'dimension')\n") ;
153 
154  if (periodicity.size()>0)
155  {
156  bool hasper=false ;
157  for (auto v : periodicity)
158  if (v)
159  hasper=true ;
160  if (!hasper)
161  periodicity.resize(0) ; // A periodic vector was defined, but with no PBC ... removing it.
162  else
163  {
164  delta.resize(dim,0) ;
165 
166  for(int i=0 ; i<dim ; i++)
167  {
168  if (periodicity[i])
169  {
170  if (boxes[i]!=1)
171  printf("WARN: using more than 1 box in the periodic dimension does not make much sense\n");
172  delta[i] = boundaries[1][i]-boundaries[0][i] ;
173  }
174  }
175  }
176  }
177 
178  if (default_density!=-1)
179  {
180  for (auto &v: files)
181  v.reader->set_default_density(default_density) ;
182  }
183  if (default_radius!=-1)
184  for (auto &v: files)
185  v.reader->set_default_radius(default_radius) ;
186 
187  for (auto &v: files)
188  v.reader->post_init() ;
189 }
190 //---------------------------------------------------
192 {
193  printf("THIS FUNCTION HAS BEEN REMOVED (identify_max_level)\n") ;
194 return -1 ;
195 }
196 //----------------------------------------------------------
197 Windows Param::identify_window(std::string windowstr)
198 {
199  if ( windowstr=="Rect3D") return Windows::Rect3D ;
200  else if ( windowstr=="Rect3DIntersect") {printf("====> DEPRECATED: misleading name, use Sphere3DIntersect instead. <=======\n") ; return Windows::Sphere3DIntersect ;}
201  else if ( windowstr=="Sphere3DIntersect") return Windows::Sphere3DIntersect ;
202  else if ( windowstr=="SphereNDIntersect") return Windows::SphereNDIntersect ;
203  else if ( windowstr=="Lucy3D") return Windows::Lucy3D ;
204  else if ( windowstr=="Hann3D") return Windows::Hann3D ;
205  else if ( windowstr=="RectND") return Windows::RectND ;
206  else if ( windowstr=="LucyND") return Windows::LucyND ;
207  else if ( windowstr=="LucyND_Periodic") return Windows::LucyND_Periodic ;
208  else if ( windowstr=="Lucy3DFancyInt") return Windows::Lucy3DFancyInt;
209  else {printf("Unknown windowing function.\n") ; return Windows::Lucy3D ; }
210 }
211 //-------------------------------------
213 {
214  for (auto & v : j)
215  {
217  if (v["type"].get<string>()=="Particle") typ = FieldType::Particle ;
218  else if (v["type"].get<string>()=="Contact") typ = FieldType::Contact ;
219  else if (v["type"].get<string>()=="Fluctuation") typ = FieldType::Fluctuation ;
220  else printf("ERR: unknown extrafields type\n") ;
221  extrafields.push_back({.name=v["name"].get<string>(), .order=static_cast<TensorOrder>(v["tensor order"].get<int>()), .type=typ}) ;
222  if (v.exist("mapping")) {extrafields.back().mapping=v["mapping"].get<string>();}
223  }
224 }
225 //---------------------------------------------------------
227 {
228  for (size_t f=0 ; f<j2.size() ; f++)
229  {
230  auto & j =j2[f] ;
231 
232  // Treat actions first if the tag is present
233  if (j.exist("action"))
234  {
235  auto v = j["action"] ;
236  if ( v.get<string>() == "remove")
237  {
238  files.erase(files.begin()+f) ;
239  continue ;
240  }
241  else if (v.get<string>() == "edit")
242  files.erase(files.begin()+f) ;
243  else if (v.get<string>() == "donothing")
244  continue ;
245  else if (v.get<string>() == "create")
246  ; // do nothing
247  else
248  printf("WARN: Unknown action command on process_file\n") ;
249  }
250 
251  FileType content = FileType::particles ;
252  if (j.exist("content"))
253  {
254  if ( j["content"].get<string>() == "particles") content = FileType::particles ;
255  else if ( j["content"].get<string>() == "contacts") content = FileType::contacts ;
256  else if ( j["content"].get<string>() == "both") content = FileType::both ;
257  else printf("WARN: unknown file content.\n") ;
258  }
259  else {printf("WARN: the key 'content' is recommended when defining a file. Set to 'particles' by default.\n") ;}
260 
261  FileFormat format ;
262  if (j.exist("format"))
263  {
264  if ( j["format"].get<string>() == "liggghts") format = FileFormat::Liggghts ;
265  else if ( j["format"].get<string>() == "NDDEM") format = FileFormat::NDDEM ;
266  else if ( j["format"].get<string>() == "mercury_legacy") format = FileFormat::MercuryData ;
267  else if ( j["format"].get<string>() == "mercury_vtu") format = FileFormat::MercuryVTU ;
268  else if ( j["format"].get<string>() == "yade") format = FileFormat::Yade ;
269  else if ( j["format"].get<string>() == "interactive") format = FileFormat::Interactive ;
270  else {printf("ERR: unknown file format.\n") ; return ;}
271  }
272  else { printf("ERR: the key 'format' is required when defining a file. \n") ; return ; }
273 
274 
275  std::map <std::string, std::string> mapping ;
276  if ( j.exist("mapping") ) mapping= j["mapping"].get<std::map <std::string, std::string>>();
277 
278  files.insert(files.begin()+f, {.fileformat=format, .filetype=content, .mapping=mapping, .reader = nullptr,}) ;
279 
280  if (files[f].fileformat == FileFormat::NDDEM)
281  files[f].reader = new NDDEMReader (j["filename"].get<string>());
282  else if (files[f].fileformat == FileFormat::Interactive)
283  files[f].reader = new InteractiveReader ();
284  #ifndef NOTALLFORMATS
285  else if (files[f].fileformat == FileFormat::Liggghts)
286  {
287  if (files[f].filetype == FileType::both)
288  files[f].reader = new LiggghtsReader (j["filename"].get<string>()) ;
289  else if (files[f].filetype == FileType::particles)
290  files[f].reader = new LiggghtsReader_particles (j["filename"].get<string>()) ;
291  else if (files[f].filetype == FileType::contacts)
292  {
293  std::vector<File>::iterator it ;
294  for (it = files.begin() ; it<files.end() ; it++)
295  if (it->filetype==FileType::particles)
296  break ;
297  if (it == files.end())
298  {
299  printf("ERR: a dump particles needs to be defined before defining a dump contacts\n") ;
300  files.pop_back() ;
301  continue ;
302  }
303  files[f].reader = new LiggghtsReader_contacts (j["filename"].get<string>(), it->reader, files[f].mapping) ;
304  }
305  }
306  else if (files[f].fileformat == FileFormat::MercuryData)
307  {
308  if (files[f].filetype == FileType::particles)
309  files[f].reader = new MercuryReader_data_particles (j["filename"].get<string>()) ;
310  else if (files[f].filetype == FileType::contacts) // TODO
311  {}//files[f].reader = new MercuryReader_contacts (files[f].path) ;
312  else
313  printf("ERR: no Mercury file format support FileType::both.\n") ;
314  }
315  else if (files[f].fileformat == FileFormat::MercuryVTU)
316  {
317  if (files[f].filetype == FileType::particles)
318  files[f].reader = new MercuryReader_vtu_particles (j["filename"].get<string>()) ;
319  //else if (files[f].filetype == FileType::contacts) // TODO
320  //{}//files[f].reader = new MercuryReader_contacts (files[f].path) ;
321  else
322  printf("ERR: no Mercury VTU file format support FileType::contacts and FileType::both.\n") ;
323  }
324  else if (files[f].fileformat == FileFormat::Yade)
325  {
326  files[f].reader = new YadeReader() ;
327  }
328  #endif
329  else
330  {
331  printf("Coarse graining has not been compiled with the requested file format") ;
332  }
333 
334  bool multifile = false ;
335  files[f].reader->path = j["filename"].get<string>() ;
336  if (j.exist("initial")) {multifile = true ; files[f].reader->filenumbering.initial = j["initial"].get<double>() ; }
337  if (j.exist("delta")) {multifile = true ; files[f].reader->filenumbering.delta = j["delta"].get<double>() ; }
338 
339  if (multifile && files[f].reader->path.find('%') == std::string::npos)
340  printf("WARN: you have provided file initial or delta numbering, but there is no pattern in your filename, which is surprising\n") ;
341  if (files[f].reader->path.find('%') != std::string::npos) files[f].reader->filenumbering.ismultifile = true ;
342 
343  }
344 }
345 //------------------
346 int Param::set_data(struct Data & D)
347 {
352 
353  D.pos.resize(dim) ; D.vel.resize(dim) ;
354  D.pospq.resize(dim) ; D.lpq.resize (dim) ;
355  D.fpq.resize (dim) ;
356  D.omega.resize(dim*(dim-1)/2);
357  D.mpq.resize (dim*(dim-1)/2) ; D.mqp.resize (dim*(dim-1)/2) ;
358 
361  for (int dd=0 ; dd<dim ; dd++)
362  {
363  D.pos[dd] = get_data(DataValue::pos, dd) ;
364  D.vel[dd] = get_data(DataValue::vel, dd) ;
365  D.pospq[dd] = get_data(DataValue::pospq, dd) ;
366  D.lpq[dd] = get_data(DataValue::lpq, dd) ;
367  D.fpq[dd] = get_data(DataValue::fpq, dd) ;
368  }
369  for (int dd=0 ; dd<dim*(dim-1)/2 ; dd++)
370  {
371  D.omega[dd] = get_data(DataValue::omega, dd) ;
372  D.mpq[dd] = get_data(DataValue::mpq, dd) ;
373  D.mqp[dd] = get_data(DataValue::mqp, dd) ;
374  }
375 
376  for (auto &v: D.extrafields)
377  {
378  for (int dd=0 ; dd<std::get<1>(v) ; dd++)
379  D.extra[std::get<2>(v)+dd] = get_data(DataValue::extra_named, dd, std::get<0>(v)) ;
380  }
381  return 0 ;
382 }
383 //------------------------------------------------------------------------
385 {
386  for (auto &v: files)
387  if (v.reader->get_num_particles() != -1)
388  return v.reader->get_num_particles() ;
389  return -1 ;
390 }
392 {
393  for (auto &v: files)
394  if (v.reader->get_num_contacts() != -1)
395  return v.reader->get_num_contacts() ;
396  return -1 ;
397 }
398 double * Param::get_data(DataValue datavalue, int dd, std::string name)
399 {
400  double * res ;
401  for (auto &v: files)
402  {
403  res = v.reader->get_data(datavalue, dd, name) ;
404  if (res != nullptr)
405  return res;
406  }
407  return nullptr ;
408 }
409 #endif
FileType
Definition: Parameters.h:8
@ particles
Definition: Parameters.h:8
@ both
Definition: Parameters.h:8
@ contacts
Definition: Parameters.h:8
FileFormat
Definition: Parameters.h:7
@ Interactive
Definition: Parameters.h:7
@ Yade
Definition: Parameters.h:7
@ MercuryVTU
Definition: Parameters.h:7
@ NDDEM
Definition: Parameters.h:7
@ Liggghts
Definition: Parameters.h:7
@ MercuryData
Definition: Parameters.h:7
Windows
Definition: WindowLibrary.h:2
@ Lucy3DFancyInt
@ LucyND_Periodic
@ SphereNDIntersect
@ Sphere3DIntersect
Definition: Reader-interactive.h:7
Definition: Reader-Liggghts.h:105
Definition: Reader-Liggghts.h:72
Definition: Reader-Liggghts.h:18
Definition: Reader-Mercury.h:56
Definition: Reader-Mercury.h:28
Definition: Reader-NDDEM.h:5
Definition: Reader.h:10
Definition: Reader-Yade.h:11
namespace for Niels Lohmann
Definition: json.hpp:20203
auto get() const noexcept(noexcept(std::declval< const basic_json_t & >().template get_impl< ValueType >(detail::priority_tag< 4 > {}))) -> decltype(std::declval< const basic_json_t & >().template get_impl< ValueType >(detail::priority_tag< 4 > {}))
get a (pointer) value (explicit)
Definition: json.hpp:21875
vector< double * > omega
Particle angular velocity.
Definition: Coarsing.h:102
AverageType
Definition: Coarsing.h:60
vector< double * > mpq
Moment of particle 1 on 2.
Definition: Coarsing.h:113
double * radius
Particle radius.
Definition: Coarsing.h:97
int N
Number of particles.
Definition: Coarsing.h:96
vector< double * > pospq
Location of contact point.
Definition: Coarsing.h:110
double * id2
Index of the second particle in contact.
Definition: Coarsing.h:109
vector< double * > fpq
Force at contact.
Definition: Coarsing.h:112
double * mass
Particle masses.
Definition: Coarsing.h:98
vector< double * > pos
Particle positions.
Definition: Coarsing.h:100
TensorOrder
Definition: Coarsing.h:58
int Ncf
Number of contacts.
Definition: Coarsing.h:107
vector< double * > extra
Definition: Coarsing.h:117
double * id1
Index of first particle in contact.
Definition: Coarsing.h:108
vector< std::tuple< std::string, int, int > > extrafields
Definition: Coarsing.h:118
FieldType
Definition: Coarsing.h:59
double * Imom
Particle moment of inertia.
Definition: Coarsing.h:99
vector< double * > lpq
Branch vector of the contact, use compute_lpq() to populate this.
Definition: Coarsing.h:111
vector< double * > mqp
Moment of particle 2 on 1.
Definition: Coarsing.h:114
vector< double * > vel
Particle velocity.
Definition: Coarsing.h:101
@ Final
Definition: Coarsing.h:60
@ None
Definition: Coarsing.h:60
@ Both
Definition: Coarsing.h:60
@ Intermediate
Definition: Coarsing.h:60
@ Particle
Definition: Coarsing.h:59
@ Fluctuation
Definition: Coarsing.h:59
@ Contact
Definition: Coarsing.h:59
DataValue
Definition: Typedefs.h:19
@ fpq
Definition: Typedefs.h:19
@ Imom
Definition: Typedefs.h:19
@ pospq
Definition: Typedefs.h:19
@ id2
Definition: Typedefs.h:19
@ lpq
Definition: Typedefs.h:19
@ pos
Definition: Typedefs.h:19
@ vel
Definition: Typedefs.h:19
@ mqp
Definition: Typedefs.h:19
@ mpq
Definition: Typedefs.h:19
@ id1
Definition: Typedefs.h:19
@ extra_named
Definition: Typedefs.h:19
@ omega
Definition: Typedefs.h:19
@ radius
Definition: Typedefs.h:19
@ mass
Definition: Typedefs.h:19
struct JsonValue json
Definition: json_parser.h:15
const GenericPointer< typename T::ValueType > T2 T::AllocatorType & a
Definition: pointer.h:1181
Data structure handling point data and contact data.
Definition: Coarsing.h:93
Definition: IOLiggghts.h:92
TensorOrder order
Definition: IOLiggghts.h:94
string name
Definition: IOLiggghts.h:93
std::optional< std::string > mapping
Definition: Parameters.h:48
int datalocation
Definition: Parameters.h:49
FieldType type
Definition: IOLiggghts.h:95
Definition: Parameters.h:34
FileFormat fileformat
Definition: Parameters.h:35
Reader * reader
Definition: Parameters.h:38
std::map< std::string, std::string > mapping
Definition: Parameters.h:37
FileType filetype
Definition: Parameters.h:36
Contains the parameters of the simulation & CG.
Definition: IOLiggghts.h:67
string windowstr
Definition: IOLiggghts.h:86
int set_data(Data &cgdata)
Definition: Parameters.h:346
int read_timestep(int ts, bool particleonly=false)
Definition: Parameters.h:116
vector< int > boxes
Definition: IOLiggghts.h:75
std::vector< File > files
Definition: Parameters.h:41
vector< vector< double > > boundaries
Definition: IOLiggghts.h:76
double default_density
Definition: Parameters.h:22
double windowsize
Definition: IOLiggghts.h:88
bool tsread
Definition: Parameters.h:54
Windows identify_window(std::string windowname)
Definition: Parameters.h:197
vector< ExtaField > extrafields
Definition: IOLiggghts.h:97
AverageType timeaverage
Definition: Parameters.h:24
vector< bool > periodicity
Definition: IOLiggghts.h:89
double * get_data(DataValue datavalue, int dd=0, std::string name="")
Definition: Parameters.h:398
vector< string > saveformat
Definition: IOLiggghts.h:78
void from_json(json &j)
string save
Definition: IOLiggghts.h:77
int maxT
Definition: IOLiggghts.h:72
void process_extrafields(json &j)
Definition: Parameters.h:212
vector< string > flags
Definition: IOLiggghts.h:73
Windows window
Definition: IOLiggghts.h:87
int get_num_particles()
Definition: Parameters.h:384
int get_num_contacts()
Definition: Parameters.h:391
void post_init()
const int dim
Definition: IOLiggghts.h:74
vector< double > delta
Definition: IOLiggghts.h:90
int identify_max_level()
Definition: Parameters.h:191
std::map< std::string, bool > requiredfieldset
Definition: Parameters.h:31
int skipT
Definition: IOLiggghts.h:71
double default_radius
Definition: Parameters.h:22
int curts
Definition: Parameters.h:53
void process_file(json &j)
Definition: Parameters.h:226