gnuplot_io.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 // C++ includes
19 #include <fstream>
20 #include <sstream>
21 #include <map>
22 
23 // Local includes
24 #include "libmesh/elem.h"
26 #include "libmesh/mesh_base.h"
27 #include "libmesh/gnuplot_io.h"
28 
29 namespace libMesh
30 {
31 
33  const std::string& title,
34  int mesh_properties)
35  :
36  MeshOutput<MeshBase> (mesh_in),
37  _title(title)
38 {
39  _grid = (mesh_properties & GRID_ON);
40  _png_output = (mesh_properties & PNG_OUTPUT);
41 }
42 
43 void GnuPlotIO::write(const std::string& fname)
44 {
45  this->write_solution(fname);
46 }
47 
48 void GnuPlotIO::write_nodal_data (const std::string& fname,
49  const std::vector<Number>& soln,
50  const std::vector<std::string>& names)
51 {
52  START_LOG("write_nodal_data()", "GnuPlotIO");
53 
54  this->write_solution(fname, &soln, &names);
55 
56  STOP_LOG("write_nodal_data()", "GnuPlotIO");
57 }
58 
59 
60 
61 
62 void GnuPlotIO::write_solution(const std::string& fname,
63  const std::vector<Number>* soln,
64  const std::vector<std::string>* names)
65 {
66  // Even when writing on a serialized ParallelMesh, we expect
67  // non-proc-0 help with calls like n_active_elem
68  // libmesh_assert_equal_to (this->mesh().processor_id(), 0);
69 
70  const MeshBase& the_mesh = MeshOutput<MeshBase>::mesh();
71 
72  dof_id_type n_active_elem = the_mesh.n_active_elem();
73 
74  if (this->mesh().processor_id() == 0)
75  {
76  std::stringstream data_stream_name;
77  data_stream_name << fname << "_data";
78  const std::string data_file_name = data_stream_name.str();
79 
80  // This class is designed only for use with 1D meshes
81  libmesh_assert_equal_to (the_mesh.mesh_dimension(), 1);
82 
83  // Make sure we have a solution to plot
84  libmesh_assert ((names != NULL) && (soln != NULL));
85 
86  // Create an output stream for script file
87  std::ofstream out_stream(fname.c_str());
88 
89  // Make sure it opened correctly
90  if (!out_stream.good())
91  libmesh_file_error(fname.c_str());
92 
93  // The number of variables in the equation system
94  const unsigned int n_vars =
95  libmesh_cast_int<unsigned int>(names->size());
96 
97  // Write header to stream
98  out_stream << "# This file was generated by gnuplot_io.C\n"
99  << "# Stores 1D solution data in GNUplot format\n"
100  << "# Execute this by loading gnuplot and typing "
101  << "\"call '" << fname << "'\"\n"
102  << "reset\n"
103  << "set title \"" << _title << "\"\n"
104  << "set xlabel \"x\"\n"
105  << "set xtics nomirror\n";
106 
107  // Loop over the elements to find the minimum and maximum x values,
108  // and also to find the element boundaries to write out as xtics
109  // if requested.
110  Real x_min=0., x_max=0.;
111 
112  // construct string for xtic positions at element edges
113  std::stringstream xtics_stream;
114 
116  const MeshBase::const_element_iterator end_it =
117  the_mesh.active_elements_end();
118 
119  unsigned int count = 0;
120 
121  for( ; it != end_it; ++it)
122  {
123  const Elem* el = *it;
124 
125  // if el is the left edge of the mesh, print its left node position
126  if(el->neighbor(0) == NULL)
127  {
128  x_min = (*(el->get_node(0)))(0);
129  xtics_stream << "\"\" " << x_min << ", \\\n";
130  }
131  if(el->neighbor(1) == NULL)
132  {
133  x_max = (*(el->get_node(1)))(0);
134  }
135  xtics_stream << "\"\" " << (*(el->get_node(1)))(0);
136 
137  if(count+1 != n_active_elem)
138  {
139  xtics_stream << ", \\\n";
140  }
141  count++;
142  }
143 
144  out_stream << "set xrange [" << x_min << ":" << x_max << "]\n";
145 
146  if(_grid)
147  out_stream << "set x2tics (" << xtics_stream.str() << ")\nset grid noxtics noytics x2tics\n";
148 
149  if(_png_output)
150  {
151  out_stream << "set terminal png\n";
152  out_stream << "set output \"" << fname << ".png\"\n";
153  }
154 
155  out_stream << "plot "
156  << axes_limits
157  << " \"" << data_file_name << "\" using 1:2 title \"" << (*names)[0]
158  << "\" with lines";
159  if(n_vars > 1)
160  {
161  for(unsigned int i=1; i<n_vars; i++)
162  {
163  out_stream << ", \\\n\"" << data_file_name << "\" using 1:" << i+2
164  << " title \"" << (*names)[i] << "\" with lines";
165  }
166  }
167 
168  out_stream.close();
169 
170 
171  // Create an output stream for data file
172  std::ofstream data(data_file_name.c_str());
173 
174  if (!data.good())
175  {
176  libMesh::err << "ERROR: opening output data file " << std::endl;
177  libmesh_error();
178  }
179 
180  // get ordered nodal data using a map
181  typedef std::pair<Real, std::vector<Number> > key_value_pair;
182  typedef std::map<Real, std::vector<Number> > map_type;
183  typedef map_type::iterator map_iterator;
184 
185  map_type node_map;
186 
187 
188  it = the_mesh.active_elements_begin();
189 
190  for ( ; it != end_it; ++it)
191  {
192  const Elem* elem = *it;
193 
194  for(unsigned int i=0; i<elem->n_nodes(); i++)
195  {
196  std::vector<Number> values;
197 
198  // Get the global id of the node
199  dof_id_type global_id = elem->node(i);
200 
201  for(unsigned int c=0; c<n_vars; c++)
202  {
203  values.push_back( (*soln)[global_id*n_vars + c] );
204  }
205 
206  node_map[ the_mesh.point(global_id)(0) ] = values;
207  }
208  }
209 
210 
211  map_iterator map_it = node_map.begin();
212  const map_iterator end_map_it = node_map.end();
213 
214  for( ; map_it != end_map_it; ++map_it)
215  {
216  key_value_pair kvp = *map_it;
217  std::vector<Number> values = kvp.second;
218 
219  data << kvp.first << "\t";
220 
221  for(unsigned int i=0; i<values.size(); i++)
222  {
223  data << values[i] << "\t";
224  }
225 
226  data << "\n";
227  }
228 
229  data.close();
230  }
231 }
232 
233 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo