MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
get-values.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 //
12 // --------------------------------------------------------------------
13 // Get Values Miniapp: Extract field values via DataCollection classes
14 // --------------------------------------------------------------------
15 //
16 // This miniapp loads previously saved data using DataCollection sub-classes,
17 // see e.g. load-dc miniapp and Example 5/5p, and output field values at a
18 // set of points. Currently, only the VisItDataCollection class is supported.
19 //
20 // Compile with: make get-values
21 //
22 // Serial sample runs:
23 // > get-values -r ../../examples/Example5 -p "0 0 0.1 0" -fn pressure
24 //
25 // Parallel sample runs:
26 // > mpirun -np 4 get-values -r ../electromagnetics/Volta-AMR-Parallel
27 // -c 1 -p "0 0 0 0.1 0.1 0.1" -fn ALL
28 //
29 // Point locations can be specified on the command line using -p or within a
30 // data file whose name can be given with option -pf. The data file format is:
31 //
32 // number_of_points space_dimension
33 // x_0 y_0 ...
34 // x_1 y_1 ...
35 // etc.
36 //
37 // By default all available fields are evaluated. The list of fields can be
38 // reduced by specifying the desired field names with -fn. The -fn option
39 // takes a space separated list of field names surrounded by quotes. Field
40 // names containing spaces, such as "Field 1" and "Field 2", can be entered as:
41 // get-values -fn "Field\ 1 Field\ 2"
42 //
43 // By default the data is written to standard out. This can be overwritten
44 // with the -o [filename] option.
45 //
46 // The output format contains comments as well as sizing information to aid in
47 // subsequent processing. The bulk of the data consists of one line per point
48 // with a 0-based integer index followed by the point coordinates and then the
49 // field data. A legend, appearing before the bulk data, shows the order of
50 // the fields along with the number of values per field (for vector data).
51 //
52 #include "mfem.hpp"
53 
54 #include <fstream>
55 #include <set>
56 #include <string>
57 
58 using namespace std;
59 using namespace mfem;
60 
61 void parseFieldNames(const char * field_name_c_str, set<string> &field_names);
62 void parsePoints(int spaceDim, const char *pts_file_c_str, Vector &pts);
63 /// @cond Suppress_Doxygen_warnings
64 void writeLegend(const DataCollection::FieldMapType &fields,
65  const set<string> & field_names, int spaceDim);
66 /// @endcond
67 
68 int main(int argc, char *argv[])
69 {
70 #ifdef MFEM_USE_MPI
71  Mpi::Init();
72  if (!Mpi::Root()) { mfem::out.Disable(); mfem::err.Disable(); }
73  Hypre::Init();
74 #endif
75 
76  // Parse command-line options.
77  const char *coll_name = NULL;
78  int cycle = 0;
79  int pad_digits_cycle = 6;
80  int pad_digits_rank = 6;
81 
82  const char *field_name_c_str = "ALL";
83  const char *pts_file_c_str = "";
84  Vector pts;
85 
86  const char *out_file_c_str = "";
87 
88  OptionsParser args(argc, argv);
89  args.AddOption(&coll_name, "-r", "--root-file",
90  "Set the VisIt data collection root file prefix.", true);
91  args.AddOption(&cycle, "-c", "--cycle", "Set the cycle index to read.");
92  args.AddOption(&pad_digits_cycle, "-pdc", "--pad-digits-cycle",
93  "Number of digits in cycle.");
94  args.AddOption(&pad_digits_rank, "-pdr", "--pad-digits-rank",
95  "Number of digits in MPI rank.");
96  args.AddOption(&pts, "-p", "--points", "List of points.");
97  args.AddOption(&field_name_c_str, "-fn", "--field-names",
98  "List of field names to get values from.");
99  args.AddOption(&pts_file_c_str, "-pf", "--point-file",
100  "Filename containing a list of points.");
101  args.AddOption(&out_file_c_str, "-o", "--output-file",
102  "Output filename.");
103  args.Parse();
104  if (!args.Good())
105  {
106  args.PrintUsage(mfem::out);
107  return 1;
108  }
109  args.PrintOptions(mfem::out);
110 
111 #ifdef MFEM_USE_MPI
112  VisItDataCollection dc(MPI_COMM_WORLD, coll_name);
113 #else
114  VisItDataCollection dc(coll_name);
115 #endif
116  dc.SetPadDigitsCycle(pad_digits_cycle);
117  dc.SetPadDigitsRank(pad_digits_rank);
118  dc.Load(cycle);
119 
120  if (dc.Error() != DataCollection::NO_ERROR)
121  {
122  mfem::out << "Error loading VisIt data collection: " << coll_name << endl;
123  return 1;
124  }
125 
126  int spaceDim = dc.GetMesh()->SpaceDimension();
127 
128  mfem::out << endl;
129  mfem::out << "Collection Name: " << dc.GetCollectionName() << endl;
130  mfem::out << "Space Dimension: " << spaceDim << endl;
131  mfem::out << "Cycle: " << dc.GetCycle() << endl;
132  mfem::out << "Time: " << dc.GetTime() << endl;
133  mfem::out << "Time Step: " << dc.GetTimeStep() << endl;
134  mfem::out << endl;
135 
136  typedef DataCollection::FieldMapType fields_t;
137  const fields_t &fields = dc.GetFieldMap();
138  // Print the names of all fields.
139  mfem::out << "fields: [ ";
140  for (fields_t::const_iterator it = fields.begin(); it != fields.end(); ++it)
141  {
142  if (it != fields.begin()) { mfem::out << ", "; }
143  mfem::out << it->first;
144  }
145  mfem::out << " ]" << endl;
146 
147  // Parsing desired field names
148  set<string> field_names;
149  parseFieldNames(field_name_c_str, field_names);
150 
151  // Print field names to be extracted
152  mfem::out << "Extracting fields: ";
153  for (set<string>::iterator it=field_names.begin();
154  it!=field_names.end(); it++)
155  {
156  mfem::out << " \"" << *it << "\"";
157  }
158  mfem::out << endl;
159 
160  parsePoints(spaceDim, pts_file_c_str, pts);
161  int npts = pts.Size() / spaceDim;
162 
163  DenseMatrix pt_mat(pts.GetData(), spaceDim, npts);
164 
165  Array<int> elem_ids;
167  int nfound = dc.GetMesh()->FindPoints(pt_mat, elem_ids, ip);
168  mfem::out << "Found " << nfound << " points." << endl;
169 
170  ofstream ofs;
171  if (strcmp(out_file_c_str,"") != 0
172 #ifdef MFEM_USE_MPI
173  && Mpi::Root()
174 #endif
175  )
176  {
177  ofs.open(out_file_c_str);
178  if (!ofs)
179  {
180  MFEM_ABORT("Failed to open output file: " << out_file_c_str << '\n');
181  }
182 
183  mfem::out.SetStream(ofs);
184  }
185 
186  // Write legend showing the order of the fields and their sizes
187  writeLegend(fields, field_names, spaceDim);
188 
189  mfem::out << "# Number of points\n" << nfound << endl;
190 
191  for (int e=0; e<elem_ids.Size(); e++)
192  {
193  if (elem_ids[e] >= 0)
194  {
195  mfem::out << e;
196  for (int d=0; d<spaceDim; d++)
197  {
198  mfem::out << ' ' << pt_mat(d, e);
199  }
200 
201  // Loop over all fields.
202  for (fields_t::const_iterator it = fields.begin();
203  it != fields.end(); ++it)
204  {
205  if (field_names.find("ALL") != field_names.end() ||
206  field_names.find(it->first) != field_names.end())
207  {
208  if (it->second->VectorDim() == 1)
209  {
210  mfem::out << ' ' << it->second->GetValue(elem_ids[e], ip[e]);
211  }
212  else
213  {
214  Vector val;
215  it->second->GetVectorValue(elem_ids[e], ip[e], val);
216  for (int d=0; d<it->second->VectorDim(); d++)
217  {
218  mfem::out << ' ' << val[d];
219  }
220  }
221  }
222  }
223  mfem::out << endl;
224  }
225  }
226  if (ofs) { ofs.close(); }
227 
228  return 0;
229 }
230 
231 void parseFieldNames(const char * field_name_c_str, set<string> &field_names)
232 {
233  string field_name_str(field_name_c_str);
234  string field_name;
235 
236  for (string::iterator it=field_name_str.begin();
237  it!=field_name_str.end(); it++)
238  {
239  if (*it == '\\')
240  {
241  it++;
242  field_name.push_back(*it);
243  }
244  else if (*it == ' ')
245  {
246  if (!field_name.empty())
247  {
248  field_names.insert(field_name);
249  }
250  field_name.clear();
251  }
252  else if (it == field_name_str.end() - 1)
253  {
254  field_name.push_back(*it);
255  field_names.insert(field_name);
256  }
257  else
258  {
259  field_name.push_back(*it);
260  }
261  }
262  if (field_names.size() == 0)
263  {
264  field_names.insert("ALL");
265  }
266 }
267 
268 void parsePoints(int spaceDim, const char *pts_file_c_str, Vector &pts)
269 {
270  if (strcmp(pts_file_c_str,"") == 0) { return; }
271 
272  ifstream ifs(pts_file_c_str);
273 
274  if (!ifs)
275  {
276  MFEM_ABORT("Failed to open point file: " << pts_file_c_str << '\n');
277  }
278 
279  int o, n, dim;
280  ifs >> n >> dim;
281 
282  MFEM_VERIFY(dim == spaceDim, "Mismatch in mesh's space dimension "
283  "and point dimension.");
284 
285  if (pts.Size() > 0 && pts.Size() % spaceDim == 0)
286  {
287  int size = pts.Size() + n * dim;
288  Vector tmp(size);
289  for (int i=0; i<pts.Size(); i++)
290  {
291  tmp[i] = pts[i];
292  }
293  o = pts.Size();
294  pts.Swap(tmp);
295  }
296  else
297  {
298  pts.SetSize(n * dim);
299  o = 0;
300  }
301 
302  for (int i=0; i<n; i++)
303  {
304  for (int d=0; d<dim; d++)
305  {
306  ifs >> pts[o + i * dim + d];
307  }
308  }
309 
310  ifs.close();
311 }
312 
314  const set<string> & field_names, int spaceDim)
315 {
316  typedef DataCollection::FieldMapType fields_t;
317 
318  // Count the number of fields to be output
319  int nfields = 1;
320  for (fields_t::const_iterator it = fields.begin();
321  it != fields.end(); ++it)
322  {
323  if (field_names.find("ALL") != field_names.end() ||
324  field_names.find(it->first) != field_names.end())
325  {
326  nfields++;
327  }
328  }
329 
330  // Write the legend showing each field name and its number of entries
331  mfem::out << "# Number of fields" << endl << nfields << endl;
332  mfem::out << "# Legend" << endl;
333  mfem::out << "# \"Index\" \"Location\":" << spaceDim;
334  for (fields_t::const_iterator it = fields.begin();
335  it != fields.end(); ++it)
336  {
337  if (field_names.find("ALL") != field_names.end() ||
338  field_names.find(it->first) != field_names.end())
339  {
340  mfem::out << " \"" << it->first << "\":" << it->second->VectorDim();
341  }
342  }
343  mfem::out << endl;
344 
345  // Write the number of entries for each field without the names
346  // which should be more convenient for parsing the output file.
347  mfem::out << spaceDim;
348  for (fields_t::const_iterator it = fields.begin();
349  it != fields.end(); ++it)
350  {
351  if (field_names.find("ALL") != field_names.end() ||
352  field_names.find(it->first) != field_names.end())
353  {
354  mfem::out << " " << it->second->VectorDim();
355  }
356  }
357  mfem::out << endl;
358 }
void writeLegend(const DataCollection::FieldMapType &fields, const set< string > &field_names, int spaceDim)
Definition: get-values.cpp:313
virtual void SetPadDigitsRank(int digits)
Set the number of digits used for the MPI rank in filenames.
double GetTime() const
Get physical time (for time-dependent simulations)
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
const std::string & GetCollectionName() const
Get the name of the collection.
void parsePoints(int spaceDim, const char *pts_file_c_str, Vector &pts)
Definition: get-values.cpp:268
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
int Error() const
Get the current error state.
void SetStream(std::ostream &os)
Replace the rdbuf() and tie() of the OutStream with that of os, enabling output.
Definition: globals.hpp:46
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
double GetTimeStep() const
Get the simulation time step (for time-dependent simulations)
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition: vector.hpp:615
const FieldMapType & GetFieldMap() const
Get a const reference to the internal field map.
int GetCycle() const
Get time cycle (for time-dependent simulations)
Data collection with VisIt I/O routines.
void parseFieldNames(const char *field_name_c_str, set< string > &field_names)
Definition: get-values.cpp:231
virtual void Load(int cycle_=0) override
Load the collection based on its VisIt data (described in its root file)
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int SpaceDimension() const
Definition: mesh.hpp:1007
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition: globals.hpp:71
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void Disable()
Disable output.
Definition: globals.hpp:54
GFieldMap::MapType FieldMapType
int dim
Definition: ex24.cpp:53
virtual void SetPadDigitsCycle(int digits)
Set the number of digits used for the cycle.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
Mesh * GetMesh()
Get a pointer to the mesh in the collection.
virtual int FindPoints(DenseMatrix &point_mat, Array< int > &elem_ids, Array< IntegrationPoint > &ips, bool warn=true, InverseElementTransformation *inv_trans=NULL)
Find the ids of the elements that contain the given points, and their corresponding reference coordin...
Definition: mesh.cpp:11898
Vector data type.
Definition: vector.hpp:60
void pts(int iphi, int t, double x[])
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
int main()
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150