1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 <ctime>
20 #include <iomanip>
21 #include <iostream>
22 #include <sstream>
23
24 // Local includes
25 #include "libmesh/postscript_io.h"
26 #include "libmesh/mesh_tools.h"
27 #include "libmesh/elem.h"
28 #include "libmesh/int_range.h"
29
30 namespace libMesh
31 {
32
33
34 // Transformation map between monomial (physical space) and Bezier bases.
35 const float PostscriptIO::_bezier_transform[3][3] =
36 {
37 {-1.f/6.f, 1.f/6.f, 1.},
38 {-1.f/6.f, 0.5, 1.f/6.f},
39 {0., 1., 0.}
40 };
41
42
PostscriptIO(const MeshBase & mesh_in)43 PostscriptIO::PostscriptIO (const MeshBase & mesh_in) :
44 MeshOutput<MeshBase> (mesh_in),
45 shade_value(0.0),
46 line_width(0.5),
47 //_M(3,3),
48 _offset(0., 0.),
49 _scale(1.0),
50 _current_point(0., 0.)
51 {
52 // This code is still undergoing some development.
53 libmesh_experimental();
54
55 // Entries of transformation matrix from physical to Bezier coords.
56 // _M(0,0) = -1./6.; _M(0,1) = 1./6.; _M(0,2) = 1.;
57 // _M(1,0) = -1./6.; _M(1,1) = 0.5 ; _M(1,2) = 1./6.;
58 // _M(2,0) = 0. ; _M(2,1) = 1. ; _M(2,2) = 0.;
59
60 // Make sure there is enough room to store Bezier coefficients.
61 _bezier_coeffs.resize(3);
62 }
63
64
65
~PostscriptIO()66 PostscriptIO::~PostscriptIO ()
67 {
68 }
69
70
71
write(const std::string & fname)72 void PostscriptIO::write (const std::string & fname)
73 {
74 // We may need to gather a DistributedMesh to output it, making that
75 // const qualifier in our constructor a dirty lie
76 MeshSerializer serialize(const_cast<MeshBase &>(this->mesh()), !_is_parallel_format);
77
78 if (this->mesh().processor_id() == 0)
79 {
80 // Get a constant reference to the mesh.
81 const MeshBase & the_mesh = MeshOutput<MeshBase>::mesh();
82
83 // Only works in 2D
84 libmesh_assert_equal_to (the_mesh.mesh_dimension(), 2);
85
86 // Create output file stream.
87 // _out is now a private member of the class.
88 _out.open(fname.c_str());
89
90 // Make sure it opened correctly
91 if (!_out.good())
92 libmesh_file_error(fname.c_str());
93
94 // The mesh bounding box gives us info about what the
95 // Postscript bounding box should be.
96 BoundingBox bbox = MeshTools::create_bounding_box(the_mesh);
97
98 // Add a little extra padding to the "true" bounding box so
99 // that we can still see the boundary
100 const Real percent_padding = 0.01;
101 const Real dx=bbox.second(0)-bbox.first(0); libmesh_assert_greater (dx, 0.0);
102 const Real dy=bbox.second(1)-bbox.first(1); libmesh_assert_greater (dy, 0.0);
103
104 const Real x_min = bbox.first(0) - percent_padding*dx;
105 const Real y_min = bbox.first(1) - percent_padding*dy;
106 const Real x_max = bbox.second(0) + percent_padding*dx;
107 const Real y_max = bbox.second(1) + percent_padding*dy;
108
109 // Width of the output as given in postscript units.
110 // This usually is given by the strange unit 1/72 inch.
111 // A width of 300 represents a size of roughly 10 cm.
112 const Real width = 300;
113 _scale = width / (x_max-x_min);
114 _offset(0) = x_min;
115 _offset(1) = y_min;
116
117 // Header writing stuff stolen from Deal.II
118 std::time_t time1= std::time (0);
119 std::tm * time = std::localtime(&time1);
120 _out << "%!PS-Adobe-2.0 EPSF-1.2" << '\n'
121 //<< "%!PS-Adobe-1.0" << '\n' // Lars' PS version
122 << "%%Filename: " << fname << '\n'
123 << "%%Title: LibMesh Output" << '\n'
124 << "%%Creator: LibMesh: A C++ finite element library" << '\n'
125 << "%%Creation Date: "
126 << time->tm_year+1900 << "/"
127 << time->tm_mon+1 << "/"
128 << time->tm_mday << " - "
129 << time->tm_hour << ":"
130 << std::setw(2) << time->tm_min << ":"
131 << std::setw(2) << time->tm_sec << '\n'
132 << "%%BoundingBox: "
133 // lower left corner
134 << "0 0 "
135 // upper right corner
136 << static_cast<unsigned int>( rint(double((x_max-x_min) * _scale )))
137 << ' '
138 << static_cast<unsigned int>( rint(double((y_max-y_min) * _scale )))
139 << '\n';
140
141 // define some abbreviations to keep
142 // the output small:
143 // m=move turtle to
144 // l=define a line
145 // s=set rgb color
146 // sg=set gray value
147 // lx=close the line and plot the line
148 // lf=close the line and fill the interior
149 _out << "/m {moveto} bind def" << '\n'
150 << "/l {lineto} bind def" << '\n'
151 << "/s {setrgbcolor} bind def" << '\n'
152 << "/sg {setgray} bind def" << '\n'
153 << "/cs {curveto stroke} bind def" << '\n'
154 << "/lx {lineto closepath stroke} bind def" << '\n'
155 << "/lf {lineto closepath fill} bind def" << '\n';
156
157 _out << "%%EndProlog" << '\n';
158 // << '\n';
159
160 // Set line width in the postscript file.
161 _out << line_width << " setlinewidth" << '\n';
162
163 // Set line cap and join options
164 _out << "1 setlinecap" << '\n';
165 _out << "1 setlinejoin" << '\n';
166
167 // allow only five digits for output (instead of the default
168 // six); this should suffice even for fine grids, but reduces
169 // the file size significantly
170 _out << std::setprecision (5);
171
172 // Loop over the active elements, draw lines for the edges. We
173 // draw even quadratic elements with straight sides, i.e. a straight
174 // line sits between each pair of vertices. Also we draw every edge
175 // for an element regardless of the fact that it may overlap with
176 // another. This would probably be a useful optimization...
177 for (const auto & elem : the_mesh.active_element_ptr_range())
178 this->plot_linear_elem(elem);
179
180 // Issue the showpage command, and we're done.
181 _out << "showpage" << std::endl;
182
183 } // end if (this->mesh().processor_id() == 0)
184 }
185
186
187
188
189
190
plot_linear_elem(const Elem * elem)191 void PostscriptIO::plot_linear_elem(const Elem * elem)
192 {
193 // Clear the string contents. Yes, this really is how you do that...
194 _cell_string.str("");
195
196 // The general strategy is:
197 // 1.) Use m := {moveto} to go to vertex 0.
198 // 2.) Use l := {lineto} commands to draw lines to vertex 1, 2, ... N-1.
199 // 3a.) Use lx := {lineto closepath stroke} command at vertex N to draw the last line.
200 // 3b.) lf := {lineto closepath fill} command to shade the cell just drawn
201 // All of our 2D elements' vertices are numbered in counterclockwise order,
202 // so we can just draw them in the same order.
203
204 // 1.)
205 _current_point = (elem->point(0) - _offset) * _scale;
206 _cell_string << _current_point(0) << " " << _current_point(1) << " "; // write x y
207 _cell_string << "m ";
208
209 // 2.)
210 const unsigned int nv=elem->n_vertices();
211 for (unsigned int v=1; v<nv-1; ++v)
212 {
213 _current_point = (elem->point(v) - _offset) * _scale;
214 _cell_string << _current_point(0) << " " << _current_point(1) << " "; // write x y
215 _cell_string << "l ";
216 }
217
218 // 3.)
219 _current_point = (elem->point(nv-1) - _offset) * _scale;
220 _cell_string << _current_point(0) << " " << _current_point(1) << " "; // write x y
221
222 // We draw the shaded (interior) parts first, if applicable.
223 if (shade_value > 0.0)
224 _out << shade_value << " sg " << _cell_string.str() << "lf\n";
225
226 // Draw the black lines (I guess we will always do this)
227 _out << "0 sg " << _cell_string.str() << "lx\n";
228 }
229
230
231
232
plot_quadratic_elem(const Elem * elem)233 void PostscriptIO::plot_quadratic_elem(const Elem * elem)
234 {
235 for (auto ns : elem->side_index_range())
236 {
237 // Build the quadratic side
238 std::unique_ptr<const Elem> side = elem->build_side_ptr(ns);
239
240 // Be sure it's quadratic (Edge2). Eventually we could
241 // handle cubic elements as well...
242 libmesh_assert_equal_to ( side->type(), EDGE3 );
243
244 _out << "0 sg ";
245
246 // Move to the first point on this side.
247 _current_point = (side->point(0) - _offset) * _scale;
248 _out << _current_point(0) << " " << _current_point(1) << " "; // write x y
249 _out << "m ";
250
251 // Compute _bezier_coeffs for this edge. This fills up
252 // the _bezier_coeffs vector.
253 this->_compute_edge_bezier_coeffs(side.get());
254
255 // Print curveto path to file
256 for (auto i : index_range(_bezier_coeffs))
257 _out << _bezier_coeffs[i](0) << " " << _bezier_coeffs[i](1) << " ";
258 _out << " cs\n";
259 }
260 }
261
262
263
264
_compute_edge_bezier_coeffs(const Elem * elem)265 void PostscriptIO::_compute_edge_bezier_coeffs(const Elem * elem)
266 {
267 // I only know how to do this for an Edge3!
268 libmesh_assert_equal_to (elem->type(), EDGE3);
269
270 // Get x-coordinates into an array, transform them,
271 // and repeat for y.
272 float phys_coords[3] = {0., 0., 0.};
273 float bez_coords[3] = {0., 0., 0.};
274
275 for (unsigned int i=0; i<2; ++i)
276 {
277 // Initialize vectors. Physical coordinates are initialized
278 // by their postscript-scaled values.
279 for (unsigned int j=0; j<3; ++j)
280 {
281 phys_coords[j] = static_cast<float>
282 ((elem->point(j)(i) - _offset(i)) * _scale);
283 bez_coords[j] = 0.; // zero out result vector
284 }
285
286 // Multiply matrix times vector
287 for (unsigned int j=0; j<3; ++j)
288 for (unsigned int k=0; k<3; ++k)
289 bez_coords[j] += _bezier_transform[j][k]*phys_coords[k];
290
291 // Store result in _bezier_coeffs
292 for (unsigned int j=0; j<3; ++j)
293 _bezier_coeffs[j](i) = phys_coords[j];
294 }
295 }
296
297 } // namespace libMesh
298