Line data    Source code

       1             : // The libMesh Finite Element Library.
       2             : // Copyright (C) 2002-2025 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"
      25             : #include "libmesh/libmesh_logging.h"
      26             : #include "libmesh/mesh_base.h"
      27             : #include "libmesh/gnuplot_io.h"
      28             : 
      29             : namespace libMesh
      30             : {
      31             : 
      32         283 : GnuPlotIO::GnuPlotIO(const MeshBase & mesh_in,
      33             :                      const std::string & title,
      34         283 :                      int mesh_properties)
      35             :   :
      36             :   MeshOutput<MeshBase> (mesh_in),
      37         283 :   _title(title)
      38             : {
      39         283 :   _grid       = (mesh_properties & GRID_ON);
      40         283 :   _png_output = (mesh_properties & PNG_OUTPUT);
      41         283 : }
      42             : 
      43           0 : void GnuPlotIO::write(const std::string & fname)
      44             : {
      45           0 :   this->write_solution(fname);
      46           0 : }
      47             : 
      48         283 : void GnuPlotIO::write_nodal_data (const std::string & fname,
      49             :                                   const std::vector<Number> & soln,
      50             :                                   const std::vector<std::string> & names)
      51             : {
      52          16 :   LOG_SCOPE("write_nodal_data()", "GnuPlotIO");
      53         283 :   this->write_solution(fname, &soln, &names);
      54         283 : }
      55             : 
      56             : 
      57             : 
      58             : 
      59         283 : void GnuPlotIO::write_solution(const std::string & fname,
      60             :                                const std::vector<Number> * soln,
      61             :                                const std::vector<std::string> * names)
      62             : {
      63             :   // Even when writing on a serialized DistributedMesh, we expect
      64             :   // non-proc-0 help with calls like n_active_elem
      65             :   // libmesh_assert_equal_to (this->mesh().processor_id(), 0);
      66             : 
      67          16 :   const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
      68             : 
      69         283 :   dof_id_type n_active_elem = the_mesh.n_active_elem();
      70             : 
      71         291 :   if (this->mesh().processor_id() == 0)
      72             :     {
      73          55 :       std::stringstream data_stream_name;
      74          47 :       data_stream_name << fname << "_data";
      75           8 :       const std::string data_file_name = data_stream_name.str();
      76             : 
      77             :       // This class is designed only for use with 1D meshes
      78           4 :       libmesh_assert_equal_to (the_mesh.mesh_dimension(), 1);
      79             : 
      80             :       // Make sure we have a solution to plot
      81           4 :       libmesh_assert ((names != nullptr) && (soln != nullptr));
      82             : 
      83             :       // Create an output stream for script file
      84          55 :       std::ofstream out_stream(fname.c_str());
      85             : 
      86             :       // Make sure it opened correctly
      87          47 :       if (!out_stream.good())
      88           0 :         libmesh_file_error(fname.c_str());
      89             : 
      90             :       // The number of variables in the equation system
      91             :       const unsigned int n_vars =
      92           8 :         cast_int<unsigned int>(names->size());
      93             : 
      94             :       // Write header to stream
      95             :       out_stream << "# This file was generated by gnuplot_io.C\n"
      96             :                  << "# Stores 1D solution data in GNUplot format\n"
      97             :                  << "# Execute this by loading gnuplot and typing "
      98             :                  << "\"call '" << fname << "'\"\n"
      99             :                  << "reset\n"
     100           4 :                  << "set title \"" << _title << "\"\n"
     101             :                  << "set xlabel \"x\"\n"
     102         125 :                  << "set xtics nomirror\n";
     103             : 
     104             :       // Loop over the elements to find the minimum and maximum x values,
     105             :       // and also to find the element boundaries to write out as xtics
     106             :       // if requested.
     107           4 :       Real x_min=0., x_max=0.;
     108             : 
     109             :       // construct string for xtic positions at element edges
     110          55 :       std::stringstream xtics_stream;
     111             : 
     112           4 :       unsigned int count = 0;
     113             : 
     114        2734 :       for (const auto & el : the_mesh.active_element_ptr_range())
     115             :         {
     116             :           // if el is the left edge of the mesh, print its left node position
     117        1570 :           if (el->neighbor_ptr(0) == nullptr)
     118             :             {
     119          51 :               x_min = (el->point(0))(0);
     120          86 :               xtics_stream << "\"\" " << x_min << ", \\\n";
     121             :             }
     122        1570 :           if (el->neighbor_ptr(1) == nullptr)
     123             :             {
     124          51 :               x_max = (el->point(1))(0);
     125             :             }
     126        1570 :           xtics_stream << "\"\" " << (el->point(1))(0);
     127             : 
     128        1446 :           if (count+1 != n_active_elem)
     129             :             {
     130        1399 :               xtics_stream << ", \\\n";
     131             :             }
     132         124 :           count++;
     133          39 :         }
     134             : 
     135         125 :       out_stream << "set xrange [" << x_min << ":" << x_max << "]\n";
     136             : 
     137          47 :       if (_grid)
     138         129 :         out_stream << "set x2tics (" << xtics_stream.str() << ")\nset grid noxtics noytics x2tics\n";
     139             : 
     140          47 :       if (_png_output)
     141             :         {
     142           0 :           out_stream << "set terminal png\n";
     143           0 :           out_stream << "set output \"" << fname << ".png\"\n";
     144             :         }
     145             : 
     146             :       out_stream << "plot "
     147           4 :                  << axes_limits
     148           8 :                  << " \"" << data_file_name << "\" using 1:2 title \"" << (*names)[0]
     149         164 :                  << "\" with lines";
     150          47 :       if (n_vars > 1)
     151             :         {
     152          24 :           for (unsigned int i=1; i<n_vars; i++)
     153             :             {
     154          22 :               out_stream << ", \\\n\"" << data_file_name << "\" using 1:" << i+2
     155          12 :                          << " title \"" << (*names)[i] << "\" with lines";
     156             :             }
     157             :         }
     158             : 
     159          47 :       out_stream.close();
     160             : 
     161             : 
     162             :       // Create an output stream for data file
     163          55 :       std::ofstream data(data_file_name.c_str());
     164             : 
     165          47 :       libmesh_error_msg_if(!data.good(), "ERROR: opening output data file " << data_file_name);
     166             : 
     167             :       // get ordered nodal data using a map
     168           8 :       std::map<Real, std::vector<Number>> node_map;
     169             : 
     170        2734 :       for (const auto & elem : the_mesh.active_element_ptr_range())
     171        5446 :         for (const auto & node : elem->node_ref_range())
     172             :           {
     173         660 :             dof_id_type global_id = node.id();
     174             : 
     175         660 :             std::vector<Number> values;
     176        8472 :             for (unsigned int c=0; c<n_vars; c++)
     177        4986 :               values.push_back((*soln)[global_id*n_vars + c]);
     178             : 
     179        3876 :             node_map[the_mesh.point(global_id)(0)] = values;
     180          39 :           }
     181             : 
     182        2524 :       for (const auto & pr : node_map)
     183             :         {
     184         210 :           const std::vector<Number> & values = pr.second;
     185             : 
     186        2477 :           data << pr.first << "\t";
     187             : 
     188        5446 :           for (const auto & val : values)
     189        2969 :             data << val << "\t";
     190             : 
     191        2477 :           data << "\n";
     192             :         }
     193             : 
     194          47 :       data.close();
     195          78 :     }
     196         283 : }
     197             : 
     198             : } // namespace libMesh