MFEM  v4.1.0
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-2020, 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 qoutes. 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 void writeLegend(const DataCollection::FieldMapType &fields,
64  const set<string> & field_names, int spaceDim);
65 
66 int main(int argc, char *argv[])
67 {
68 #ifdef MFEM_USE_MPI
69  MPI_Session mpi;
70  if (!mpi.Root()) { mfem::out.Disable(); mfem::err.Disable(); }
71 #endif
72 
73  // Parse command-line options.
74  const char *coll_name = NULL;
75  int cycle = 0;
76 
77  const char *field_name_c_str = "ALL";
78  const char *pts_file_c_str = "";
79  Vector pts;
80 
81  const char *out_file_c_str = "";
82 
83  OptionsParser args(argc, argv);
84  args.AddOption(&coll_name, "-r", "--root-file",
85  "Set the VisIt data collection root file prefix.", true);
86  args.AddOption(&cycle, "-c", "--cycle", "Set the cycle index to read.");
87  args.AddOption(&pts, "-p", "--points", "List of points.");
88  args.AddOption(&field_name_c_str, "-fn", "--field-names",
89  "List of field names to get values from.");
90  args.AddOption(&pts_file_c_str, "-pf", "--point-file",
91  "Filename containing a list of points.");
92  args.AddOption(&out_file_c_str, "-o", "--output-file",
93  "Output filename.");
94  args.Parse();
95  if (!args.Good())
96  {
97  args.PrintUsage(mfem::out);
98  return 1;
99  }
100  args.PrintOptions(mfem::out);
101 
102 #ifdef MFEM_USE_MPI
103  VisItDataCollection dc(MPI_COMM_WORLD, coll_name);
104 #else
105  VisItDataCollection dc(coll_name);
106 #endif
107  dc.Load(cycle);
108 
109  if (dc.Error() != DataCollection::NO_ERROR)
110  {
111  mfem::out << "Error loading VisIt data collection: " << coll_name << endl;
112  return 1;
113  }
114 
115  int spaceDim = dc.GetMesh()->SpaceDimension();
116 
117  mfem::out << endl;
118  mfem::out << "Collection Name: " << dc.GetCollectionName() << endl;
119  mfem::out << "Space Dimension: " << spaceDim << endl;
120  mfem::out << "Cycle: " << dc.GetCycle() << endl;
121  mfem::out << "Time: " << dc.GetTime() << endl;
122  mfem::out << "Time Step: " << dc.GetTimeStep() << endl;
123  mfem::out << endl;
124 
125  typedef DataCollection::FieldMapType fields_t;
126  const fields_t &fields = dc.GetFieldMap();
127  // Print the names of all fields.
128  mfem::out << "fields: [ ";
129  for (fields_t::const_iterator it = fields.begin(); it != fields.end(); ++it)
130  {
131  if (it != fields.begin()) { mfem::out << ", "; }
132  mfem::out << it->first;
133  }
134  mfem::out << " ]" << endl;
135 
136  // Parsing desired field names
137  set<string> field_names;
138  parseFieldNames(field_name_c_str, field_names);
139 
140  // Print field names to be extracted
141  mfem::out << "Extracting fields: ";
142  for (set<string>::iterator it=field_names.begin();
143  it!=field_names.end(); it++)
144  {
145  mfem::out << " \"" << *it << "\"";
146  }
147  mfem::out << endl;
148 
149  parsePoints(spaceDim, pts_file_c_str, pts);
150  int npts = pts.Size() / spaceDim;
151 
152  DenseMatrix pt_mat(pts.GetData(), spaceDim, npts);
153 
154  Array<int> elem_ids;
156  int nfound = dc.GetMesh()->FindPoints(pt_mat, elem_ids, ip);
157  mfem::out << "Found " << nfound << " points." << endl;
158 
159  ofstream ofs;
160  if (strcmp(out_file_c_str,"") != 0
161 #ifdef MFEM_USE_MPI
162  && mpi.Root()
163 #endif
164  )
165  {
166  ofs.open(out_file_c_str);
167  if (!ofs)
168  {
169  MFEM_ABORT("Failed to open output file: " << out_file_c_str << '\n');
170  }
171 
172  mfem::out.SetStream(ofs);
173  }
174 
175  // Write legend showing the order of the fields and their sizes
176  writeLegend(fields, field_names, spaceDim);
177 
178  mfem::out << "# Number of points\n" << nfound << endl;
179 
180  for (int e=0; e<elem_ids.Size(); e++)
181  {
182  if (elem_ids[e] >= 0)
183  {
184  mfem::out << e;
185  for (int d=0; d<spaceDim; d++)
186  {
187  mfem::out << ' ' << pt_mat(d, e);
188  }
189 
190  // Loop over all fields.
191  for (fields_t::const_iterator it = fields.begin();
192  it != fields.end(); ++it)
193  {
194  if (field_names.find("ALL") != field_names.end() ||
195  field_names.find(it->first) != field_names.end())
196  {
197  if (it->second->VectorDim() == 1)
198  {
199  mfem::out << ' ' << it->second->GetValue(elem_ids[e], ip[e]);
200  }
201  else
202  {
203  Vector val;
204  it->second->GetVectorValue(elem_ids[e], ip[e], val);
205  for (int d=0; d<it->second->VectorDim(); d++)
206  {
207  mfem::out << ' ' << val[d];
208  }
209  }
210  }
211  }
212  mfem::out << endl;
213  }
214  }
215  if (ofs) { ofs.close(); }
216 
217  return 0;
218 }
219 
220 void parseFieldNames(const char * field_name_c_str, set<string> &field_names)
221 {
222  string field_name_str(field_name_c_str);
223  string field_name;
224 
225  for (string::iterator it=field_name_str.begin();
226  it!=field_name_str.end(); it++)
227  {
228  if (*it == '\\')
229  {
230  it++;
231  field_name.push_back(*it);
232  }
233  else if (*it == ' ')
234  {
235  if (!field_name.empty())
236  {
237  field_names.insert(field_name);
238  }
239  field_name.clear();
240  }
241  else if (it == field_name_str.end() - 1)
242  {
243  field_name.push_back(*it);
244  field_names.insert(field_name);
245  }
246  else
247  {
248  field_name.push_back(*it);
249  }
250  }
251  if (field_names.size() == 0)
252  {
253  field_names.insert("ALL");
254  }
255 }
256 
257 void parsePoints(int spaceDim, const char *pts_file_c_str, Vector &pts)
258 {
259  if (strcmp(pts_file_c_str,"") == 0) { return; }
260 
261  ifstream ifs(pts_file_c_str);
262 
263  if (!ifs)
264  {
265  MFEM_ABORT("Failed to open point file: " << pts_file_c_str << '\n');
266  }
267 
268  int o, n, dim;
269  ifs >> n >> dim;
270 
271  MFEM_VERIFY(dim == spaceDim, "Mismatch in mesh's space dimension "
272  "and point dimension.");
273 
274  if (pts.Size() > 0 && pts.Size() % spaceDim == 0)
275  {
276  int size = pts.Size() + n * dim;
277  Vector tmp(size);
278  for (int i=0; i<pts.Size(); i++)
279  {
280  tmp[i] = pts[i];
281  }
282  o = pts.Size();
283  pts.Swap(tmp);
284  }
285  else
286  {
287  pts.SetSize(n * dim);
288  o = 0;
289  }
290 
291  for (int i=0; i<n; i++)
292  {
293  for (int d=0; d<dim; d++)
294  {
295  ifs >> pts[o + i * dim + d];
296  }
297  }
298 
299  ifs.close();
300 }
301 
303  const set<string> & field_names, int spaceDim)
304 {
305  typedef DataCollection::FieldMapType fields_t;
306 
307  // Count the number of fields to be output
308  int nfields = 1;
309  for (fields_t::const_iterator it = fields.begin();
310  it != fields.end(); ++it)
311  {
312  if (field_names.find("ALL") != field_names.end() ||
313  field_names.find(it->first) != field_names.end())
314  {
315  nfields++;
316  }
317  }
318 
319  // Write the legend showing each field name and its number of entries
320  mfem::out << "# Number of fields" << endl << nfields << endl;
321  mfem::out << "# Legend" << endl;
322  mfem::out << "# \"Index\" \"Location\":" << spaceDim;
323  for (fields_t::const_iterator it = fields.begin();
324  it != fields.end(); ++it)
325  {
326  if (field_names.find("ALL") != field_names.end() ||
327  field_names.find(it->first) != field_names.end())
328  {
329  mfem::out << " \"" << it->first << "\":" << it->second->VectorDim();
330  }
331  }
332  mfem::out << endl;
333 
334  // Write the number of entries for each field without the names
335  // which should be more convenient for parsing the output file.
336  mfem::out << spaceDim;
337  for (fields_t::const_iterator it = fields.begin();
338  it != fields.end(); ++it)
339  {
340  if (field_names.find("ALL") != field_names.end() ||
341  field_names.find(it->first) != field_names.end())
342  {
343  mfem::out << " " << it->second->VectorDim();
344  }
345  }
346  mfem::out << endl;
347 }
void writeLegend(const DataCollection::FieldMapType &fields, const set< string > &field_names, int spaceDim)
Definition: get-values.cpp:302
double GetTime() const
Get physical time (for time-dependent simulations)
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
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:257
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
int Error() const
Get the current error state.
int main(int argc, char *argv[])
Definition: ex1.cpp:62
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:166
double GetTimeStep() const
Get the simulation time step (for time-dependent simulations)
A simple convenience class that calls MPI_Init() at construction and MPI_Finalize() at destruction...
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition: vector.hpp:503
void SetStream(std::ostream &out)
Replace the rdbuf() and tie() of the OutStream with that of out, enabling output. ...
Definition: globals.hpp:46
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:220
bool Root() const
Return true if WorldRank() == 0.
void PrintUsage(std::ostream &out) const
Definition: optparser.cpp:434
virtual void Load(int cycle_=0)
Load the collection based on its VisIt data (described in its root file)
int SpaceDimension() const
Definition: mesh.hpp:745
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)
Definition: optparser.hpp:76
void Disable()
Disable output.
Definition: globals.hpp:54
GFieldMap::MapType FieldMapType
int dim
Definition: ex24.cpp:43
void PrintOptions(std::ostream &out) const
Definition: optparser.cpp:304
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:10133
Vector data type.
Definition: vector.hpp:48
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
bool Good() const
Definition: optparser.hpp:125