Branch data Line data Source code
1 : : /* *****************************************************************
2 : : MESQUITE -- The Mesh Quality Improvement Toolkit
3 : :
4 : : Copyright 2005 Lawrence Livermore National Laboratory. Under
5 : : the terms of Contract B545069 with the University of Wisconsin --
6 : : Madison, Lawrence Livermore National Laboratory retains certain
7 : : rights in this software.
8 : :
9 : : This library is free software; you can redistribute it and/or
10 : : modify it under the terms of the GNU Lesser General Public
11 : : License as published by the Free Software Foundation; either
12 : : version 2.1 of the License, or (at your option) any later version.
13 : :
14 : : This library is distributed in the hope that it will be useful,
15 : : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 : : Lesser General Public License for more details.
18 : :
19 : : You should have received a copy of the GNU Lesser General Public License
20 : : (lgpl.txt) along with this library; if not, write to the Free Software
21 : : Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
22 : :
23 : : [email protected]
24 : :
25 : : ***************************************************************** */
26 : :
27 : : #ifndef MSQ_MESH_WRITER_CPP
28 : : #define MSQ_MESH_WRITER_CPP
29 : :
30 : : #include "MeshWriter.hpp"
31 : : #include "Mesquite.hpp"
32 : : #include "MeshImpl.hpp"
33 : : #include "MsqError.hpp"
34 : : #include "PatchData.hpp"
35 : : #include "PlanarDomain.hpp"
36 : : #include "VtkTypeInfo.hpp"
37 : : #include "EdgeIterator.hpp"
38 : :
39 : : #include <memory>
40 : : #include <limits>
41 : : #include <vector>
42 : : #include <algorithm>
43 : : #include <fstream>
44 : : #include <string>
45 : : #include <iomanip>
46 : : using namespace std;
47 : :
48 : : #include <stdio.h>
49 : :
50 : : namespace MBMesquite
51 : : {
52 : :
53 : : namespace MeshWriter
54 : : {
55 : :
56 : : /**\brief Transform from coordinates in the XY-plane
57 : : * to graphics coordinates.
58 : : */
59 : : class Transform2D
60 : : {
61 : : public:
62 : : Transform2D( PatchData* pd, Projection& proj, unsigned width, unsigned height, bool flip_about_horizontal );
63 : :
64 : : Transform2D( const Vector3D* verts, size_t num_vert, Projection& projection, unsigned width, unsigned height );
65 : :
66 : : void transform( const Vector3D& coords, int& horizontal, int& vertical ) const;
67 : :
68 : 0 : int max_horizontal() const
69 : : {
70 : 0 : return horizMax;
71 : : }
72 : 0 : int max_vertical() const
73 : : {
74 : 0 : return vertMax;
75 : : }
76 : :
77 : : private:
78 : : Projection& myProj;
79 : : float myScale;
80 : : int horizOffset, vertOffset;
81 : : int horizMax, vertMax;
82 : : int vertSign;
83 : : };
84 : :
85 : : /* Write VTK file
86 : : *
87 : : * Copied from src/Mesh/MeshSet.cpp and adapted for removal of
88 : : * MeshSet class by J.Kraftcheck on 2005-7-28.
89 : : *
90 : : * This code is provided mainly for debugging. A more efficient
91 : : * and complete writer implementation is provided in the MeshImpl
92 : : * class for saving meshes that were read from a file initially.
93 : : */
94 : 3 : void write_vtk( Mesh* mesh, const char* out_filename, MsqError& err )
95 : : {
96 [ - + ][ - + ]: 3 : if( MeshImpl* msq_mesh = dynamic_cast< MeshImpl* >( mesh ) )
97 : : {
98 [ # # ][ # # ]: 0 : msq_mesh->write_vtk( out_filename, err );MSQ_CHKERR( err );
[ # # ][ # # ]
[ # # ]
99 : 3 : return;
100 : : }
101 : :
102 : : // loads a global patch
103 [ + - ]: 3 : PatchData pd;
104 [ + - ]: 3 : pd.set_mesh( mesh );
105 [ + - ][ + - ]: 3 : pd.fill_global_patch( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
106 : :
107 : : // write mesh
108 [ + - ][ + - ]: 3 : write_vtk( pd, out_filename, err );MSQ_CHKERR( err );
[ - + ][ # # ]
[ # # ][ + - ]
109 : : }
110 : :
111 : 3 : void write_vtk( PatchData& pd, const char* out_filename, MsqError& err, const Vector3D* OF_gradient )
112 : : {
113 : : // Open the file
114 [ + - ][ + - ]: 3 : std::ofstream file( out_filename );
115 [ + - ][ - + ]: 3 : if( !file )
116 : : {
117 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
118 : 0 : return;
119 : : }
120 : :
121 : : // Write a header
122 [ + - ]: 3 : file << "# vtk DataFile Version 2.0\n";
123 [ + - ][ + - ]: 3 : file << "Mesquite Mesh " << out_filename << " .\n";
[ + - ]
124 [ + - ]: 3 : file << "ASCII\n";
125 [ + - ]: 3 : file << "DATASET UNSTRUCTURED_GRID\n";
126 : :
127 : : // Write vertex coordinates
128 [ + - ][ + - ]: 3 : file << "POINTS " << pd.num_nodes() << " float\n";
[ + - ][ + - ]
129 : : size_t i;
130 [ + - ][ + + ]: 3996 : for( i = 0; i < pd.num_nodes(); i++ )
131 : : {
132 [ + - ][ + - ]: 3993 : file << pd.vertex_by_index( i )[0] << ' ' << pd.vertex_by_index( i )[1] << ' ' << pd.vertex_by_index( i )[2]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
133 [ + - ]: 3993 : << '\n';
134 : : }
135 : :
136 : : // Write out the connectivity table
137 : 3 : size_t connectivity_size = 0;
138 [ + - ][ + + ]: 3003 : for( i = 0; i < pd.num_elements(); ++i )
139 [ + - ][ + - ]: 3000 : connectivity_size += pd.element_by_index( i ).node_count() + 1;
140 : :
141 [ + - ][ + - ]: 3 : file << "CELLS " << pd.num_elements() << ' ' << connectivity_size << '\n';
[ + - ][ + - ]
[ + - ][ + - ]
142 [ + - ][ + + ]: 3003 : for( i = 0; i < pd.num_elements(); i++ )
143 : : {
144 [ + - ]: 3000 : std::vector< size_t > vtx_indices;
145 [ + - ][ + - ]: 3000 : pd.element_by_index( i ).get_node_indices( vtx_indices );
146 : :
147 : : // Convert native to VTK node order, if not the same
148 : : const VtkTypeInfo* info =
149 [ + - ][ + - ]: 3000 : VtkTypeInfo::find_type( pd.element_by_index( i ).get_element_type(), vtx_indices.size(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ - + ][ # # ]
[ # # ][ - + ]
150 [ + - ]: 3000 : info->mesquiteToVtkOrder( vtx_indices );
151 : :
152 [ + - ]: 3000 : file << vtx_indices.size();
153 [ + + ]: 27000 : for( std::size_t j = 0; j < vtx_indices.size(); ++j )
154 : : {
155 [ + - ][ + - ]: 24000 : file << ' ' << vtx_indices[j];
[ + - ]
156 : : }
157 [ + - ][ + - ]: 3000 : file << '\n';
158 : 3000 : }
159 : :
160 : : // Write out the element types
161 [ + - ][ + - ]: 3 : file << "CELL_TYPES " << pd.num_elements() << '\n';
[ + - ][ + - ]
162 [ + - ][ + + ]: 3003 : for( i = 0; i < pd.num_elements(); i++ )
163 : : {
164 [ + - ]: 3000 : const VtkTypeInfo* info = VtkTypeInfo::find_type( pd.element_by_index( i ).get_element_type(),
165 [ + - ][ + - ]: 6000 : pd.element_by_index( i ).node_count(), err );MSQ_ERRRTN( err );
[ + - ][ + - ]
[ + - ][ - + ]
[ # # ][ # # ]
[ - + ]
166 [ + - ][ + - ]: 3000 : file << info->vtkType << '\n';
167 : : }
168 : :
169 : : // Write out which points are fixed.
170 [ + - ][ + - ]: 3 : file << "POINT_DATA " << pd.num_nodes() << "\nSCALARS fixed int\nLOOKUP_TABLE default\n";
[ + - ][ + - ]
171 [ + - ][ + + ]: 3996 : for( i = 0; i < pd.num_nodes(); ++i )
172 : : {
173 [ + - ][ + - ]: 3993 : if( pd.vertex_by_index( i ).get_flags() & MsqVertex::MSQ_CULLED )
[ - + ]
174 [ # # ]: 0 : file << "1\n";
175 : : else
176 [ + - ]: 3993 : file << "0\n";
177 : : }
178 [ + - ]: 3 : file << "SCALARS culled short\nLOOKUP_TABLE default\n";
179 [ + - ][ + + ]: 3996 : for( i = 0; i < pd.num_nodes(); ++i )
180 : : {
181 [ + - ][ + - ]: 3993 : if( pd.vertex_by_index( i ).is_free_vertex() )
[ + + ]
182 [ + - ]: 2789 : file << "0\n";
183 : : else
184 [ + - ]: 1204 : file << "1\n";
185 : : }
186 : :
187 [ - + ]: 3 : if( OF_gradient )
188 : : {
189 [ # # ]: 0 : file << "VECTORS gradient double\n";
190 [ # # ][ # # ]: 0 : for( i = 0; i < pd.num_free_vertices(); ++i )
191 [ # # ][ # # ]: 0 : file << OF_gradient[i].x() << " " << OF_gradient[i].y() << " " << OF_gradient[i].z() << "\n";
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
192 [ # # ][ # # ]: 0 : for( i = pd.num_free_vertices(); i < pd.num_nodes(); ++i )
[ # # ]
193 [ # # ]: 0 : file << "0.0 0.0 0.0\n";
194 : : }
195 : :
196 : : // Close the file
197 [ + - ][ + - ]: 3 : file.close();
198 : : }
199 : :
200 : 11 : void write_gnuplot( Mesh* mesh, const char* out_filebase, MsqError& err )
201 : : {
202 : : // loads a global patch
203 [ + - ]: 11 : PatchData pd;
204 [ + - ]: 11 : pd.set_mesh( mesh );
205 [ + - ][ + - ]: 22 : pd.fill_global_patch( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
206 [ + - ][ + - ]: 11 : write_gnuplot( pd, out_filebase, err );
207 : : }
208 : :
209 : 0 : void write_gnuplot( Mesh* mesh, std::vector< Mesh::ElementHandle >& elems, const char* out_filebase, MsqError& err )
210 : : {
211 : : // loads a global patch
212 [ # # ]: 0 : PatchData pd;
213 [ # # ]: 0 : pd.set_mesh( mesh );
214 [ # # ][ # # ]: 0 : std::vector< Mesh::VertexHandle > verts;
215 [ # # ][ # # ]: 0 : pd.set_mesh_entities( elems, verts, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
216 [ # # ][ # # ]: 0 : write_gnuplot( pd, out_filebase, err );
217 : : }
218 : :
219 : : /* Writes a gnuplot file directly from the MeshSet.
220 : : * This means that any mesh imported successfully into Mesquite
221 : : * can be outputed in gnuplot format.
222 : : *
223 : : * Within gnuplot, use \b plot 'file1.gpt' w l, 'file2.gpt' w l
224 : : *
225 : : * This is not geared for performance, since it has to load a global Patch from
226 : : * the mesh to write a mesh file.
227 : : *
228 : : * Copied from src/Mesh/MeshSet.cpp and adapted for removal of
229 : : * MeshSet class by J.Kraftcheck on 2005-7-28.
230 : : *
231 : : * Re-written to use EdgeIterator by J.Kraftcheck on 2009-6-11
232 : : */
233 : 11 : void write_gnuplot( PatchData& pd, const char* out_filebase, MsqError& err )
234 : : {
235 : : // Open the file
236 [ + - ]: 11 : std::string out_filename = out_filebase;
237 [ + - ]: 11 : out_filename += ".gpt";
238 [ + - ][ + - ]: 22 : std::ofstream file( out_filename.c_str() );
[ + - ]
239 [ + - ][ - + ]: 11 : if( !file )
240 : : {
241 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
242 : 0 : return;
243 : : }
244 : :
245 [ + - ][ + - ]: 22 : EdgeIterator edges( &pd, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
[ + - ]
246 : :
247 : : // Write a header
248 [ + - ]: 11 : file << "\n";
249 : :
250 [ + - ][ + + ]: 143 : while( !edges.is_at_end() )
251 : : {
252 [ + - ]: 132 : const Vector3D& s = edges.start();
253 [ + - ]: 132 : const Vector3D& e = edges.end();
254 [ + - ]: 132 : const Vector3D* m = edges.mid();
255 : :
256 [ + - ][ + - ]: 132 : file << s[0] << ' ' << s[1] << ' ' << s[2] << std::endl;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
257 [ - + ][ # # ]: 132 : if( m ) file << ( *m )[0] << ' ' << ( *m )[1] << ' ' << ( *m )[2] << std::endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
258 [ + - ][ + - ]: 132 : file << e[0] << ' ' << e[1] << ' ' << e[2] << std::endl;
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ][ + - ]
[ + - ]
259 [ + - ][ + - ]: 132 : file << std::endl << std::endl;
260 : :
261 [ + - ][ + - ]: 132 : edges.step( err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
262 : : }
263 : :
264 : : // Close the file
265 [ + - ][ + - ]: 22 : file.close();
266 : : }
267 : :
268 : : /** Helper function for write_gnuplot_animator and write_gnuplot_overlay
269 : : *
270 : : * Read a set of input files to determine the bounding box
271 : : * of the combined data.
272 : : */
273 : 1 : static void find_gnuplot_agregate_range( int count, const char* basename, Vector3D& min, Vector3D& max,
274 : : MsqError& err )
275 : : {
276 : : // read all input files to determine extents
277 [ + - ]: 1 : min = Vector3D( HUGE_VAL, HUGE_VAL, HUGE_VAL );
278 [ + - ]: 1 : max = Vector3D( -HUGE_VAL, -HUGE_VAL, -HUGE_VAL );
279 [ + + ]: 12 : for( int i = 0; i <= count; ++i )
280 : : {
281 [ + - ][ + - ]: 11 : stringstream s;
282 [ + - ][ + - ]: 11 : s << basename << '.' << i << ".gpt";
[ + - ][ + - ]
283 [ + - ][ + - ]: 22 : ifstream infile( s.str().c_str() );
[ + - ]
284 [ + - ][ - + ]: 11 : if( !infile )
285 : : {
286 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( s.str(), MsqError::FILE_ACCESS );
[ # # ]
287 : 1 : return;
288 : : }
289 : : double c[3];
290 [ + - ][ + - ]: 275 : while( infile >> c[0] >> c[1] >> c[2] )
[ + - ][ + - ]
[ + + ][ + - ]
291 : : {
292 [ + + ]: 1056 : for( int j = 0; j < 3; ++j )
293 : : {
294 [ + - ][ + + ]: 792 : if( c[j] < min[j] ) min[j] = c[j];
[ + - ]
295 [ + - ][ + + ]: 792 : if( c[j] > max[j] ) max[j] = c[j];
[ + - ]
296 : : }
297 : : }
298 : 11 : }
299 : : }
300 : :
301 : : /** Write a GNU Plot script to produce an animation from a
302 : : * sequence of data files
303 : : */
304 : 0 : void write_gnuplot_animator( int count, const char* basename, MsqError& err )
305 : : {
306 [ # # ]: 0 : if( count <= 0 ) return;
307 : :
308 : 0 : const int DELAY = 10;
309 : 0 : const int WIDTH = 640;
310 : 0 : const int HEIGHT = 480;
311 : :
312 : : // read all input files to determine extents
313 [ # # ][ # # ]: 0 : Vector3D min, max;
314 [ # # ][ # # ]: 0 : find_gnuplot_agregate_range( count, basename, min, max, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
315 : :
316 : : // chose coordinate plane to plot in
317 [ # # ]: 0 : Vector3D range = max - min;
318 : 0 : int haxis = 0, vaxis = 1;
319 [ # # ][ # # ]: 0 : if( range[0] < range[1] && range[1] < range[2] )
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
320 : : {
321 : 0 : haxis = 1;
322 : 0 : vaxis = 2;
323 : : }
324 [ # # ][ # # ]: 0 : else if( range[1] < range[2] )
[ # # ]
325 : : {
326 : 0 : vaxis = 2;
327 : : }
328 : :
329 : : // open output file
330 [ # # ]: 0 : string base( basename );
331 [ # # ][ # # ]: 0 : ofstream file( ( string( basename ) + ".gnuplot" ).c_str() );
[ # # ][ # # ]
[ # # ]
332 [ # # ][ # # ]: 0 : if( !file )
333 : : {
334 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
335 : 0 : return;
336 : : }
337 : :
338 : : // write header
339 [ # # ][ # # ]: 0 : file << "#!gnuplot" << endl;
340 [ # # ][ # # ]: 0 : file << "#" << endl;
341 [ # # ][ # # ]: 0 : file << "# Mesquite Animation of " << basename << ".0 to " << basename << '.' << count << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
342 [ # # ][ # # ]: 0 : file << "#" << endl;
343 : :
344 : : // write plot settings
345 [ # # ][ # # ]: 0 : file << "set term gif animate transparent opt delay " << DELAY << " size " << WIDTH << "," << HEIGHT << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
346 [ # # ][ # # ]: 0 : file << "set xrange [" << min[haxis] - 0.05 * range[haxis] << ":" << max[haxis] + 0.05 * range[haxis] << "]"
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
347 [ # # ]: 0 : << endl;
348 [ # # ][ # # ]: 0 : file << "set yrange [" << min[vaxis] - 0.05 * range[vaxis] << ":" << max[vaxis] + 0.05 * range[vaxis] << "]"
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
349 [ # # ]: 0 : << endl;
350 [ # # ][ # # ]: 0 : file << "set title '" << basename << "'" << endl;
[ # # ][ # # ]
351 [ # # ][ # # ]: 0 : file << "set output '" << basename << ".gif'" << endl;
[ # # ][ # # ]
352 : :
353 : : // plot data
354 [ # # ][ # # ]: 0 : for( int i = 0; i <= count; ++i )
355 [ # # ][ # # ]: 0 : file << "plot '" << basename << '.' << i << ".gpt'"
[ # # ][ # # ]
[ # # ]
356 [ # # ][ # # ]: 0 : << " using " << haxis + 1 << ":" << vaxis + 1 << " w l" << endl;
[ # # ][ # # ]
[ # # ][ # # ]
357 : : }
358 : :
359 : 11 : static unsigned red( int i, int c )
360 : : {
361 [ + + ]: 11 : if( i * 4 <= c )
362 : 3 : return 255;
363 [ + + ]: 8 : else if( i * 4 >= 3 * c )
364 : 3 : return 0;
365 : : else
366 : 5 : return 384 - i * 511 / c;
367 : : }
368 : :
369 : 11 : static unsigned green( int i, int c )
370 : : {
371 [ + + ]: 11 : if( i * 4 < c )
372 : 3 : return i * 510 / c;
373 [ + + ]: 8 : else if( i * 4 > 3 * c )
374 : 3 : return 1023 - i * 1023 / c;
375 : : else
376 : 5 : return 255;
377 : : }
378 : :
379 : 11 : static unsigned blue( int i, int c )
380 : : {
381 [ + + ]: 11 : if( i * 4 <= c )
382 : 3 : return 0;
383 [ + + ]: 8 : else if( i * 4 >= 3 * c )
384 : 3 : return 255;
385 : : else
386 : 5 : return i * 511 / c - 127;
387 : : }
388 : :
389 : : /** Write a GNU Plot script to produce a single plot from a
390 : : * sequence of data files
391 : : */
392 : 1 : void write_gnuplot_overlay( int count, const char* basename, MsqError& err )
393 : : {
394 [ - + ]: 2 : if( count <= 0 ) return;
395 : :
396 : 1 : const int WIDTH = 640;
397 : 1 : const int HEIGHT = 480;
398 : :
399 : : // read all input files to determine extents
400 [ + - ][ + - ]: 1 : Vector3D min, max;
401 [ + - ][ + - ]: 1 : find_gnuplot_agregate_range( count, basename, min, max, err );MSQ_ERRRTN( err );
[ - + ][ # # ]
[ # # ][ - + ]
402 : :
403 : : // chose coordinate plane to plot in
404 [ + - ]: 1 : Vector3D range = max - min;
405 : 1 : int haxis = 0, vaxis = 1;
406 [ + - ][ + - ]: 1 : if( range[0] < range[1] && range[1] < range[2] )
[ - + ][ # # ]
[ # # ][ # # ]
[ - + ]
407 : : {
408 : 0 : haxis = 1;
409 : 0 : vaxis = 2;
410 : : }
411 [ + - ][ + - ]: 1 : else if( range[1] < range[2] )
[ - + ]
412 : : {
413 : 0 : vaxis = 2;
414 : : }
415 : :
416 : : // open output file
417 [ + - ]: 1 : string base( basename );
418 [ + - ][ + - ]: 1 : FILE* file = fopen( ( string( basename ) + ".gnuplot" ).c_str(), "w" );
[ + - ]
419 [ - + ]: 1 : if( !file )
420 : : {
421 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
422 : 0 : return;
423 : : }
424 : :
425 : : // write header
426 [ + - ]: 1 : fprintf( file, "#!gnuplot\n" );
427 [ + - ]: 1 : fprintf( file, "#\n" );
428 [ + - ]: 1 : fprintf( file, "# Mesquite Overlay of %s.0 to %s.%d\n", basename, basename, count );
429 [ + - ]: 1 : fprintf( file, "#\n" );
430 : :
431 : : // write plot settings
432 [ + - ]: 1 : fprintf( file, "set term gif size %d,%d\n", WIDTH, HEIGHT );
433 [ + - ][ + - ]: 1 : fprintf( file, "set xrange [%f:%f]\n", min[haxis] - 0.05 * range[haxis], max[haxis] + 0.05 * range[haxis] );
[ + - ][ + - ]
[ + - ]
434 [ + - ][ + - ]: 1 : fprintf( file, "set yrange [%f:%f]\n", min[vaxis] - 0.05 * range[vaxis], max[vaxis] + 0.05 * range[vaxis] );
[ + - ][ + - ]
[ + - ]
435 [ + - ]: 1 : fprintf( file, "set title '%s'\n", basename );
436 [ + - ]: 1 : fprintf( file, "set output '%s.gif'\n", basename );
437 : :
438 : : // plot data
439 : : fprintf( file, "plot '%s.0.gpt' using %d:%d w l lc rgb '#%02x%02x%02x' title 't0'", basename, haxis + 1,
440 [ + - ]: 1 : vaxis + 1, red( 0, count ), green( 0, count ), blue( 0, count ) );
441 [ + + ]: 11 : for( int i = 1; i <= count; ++i )
442 : : {
443 : : fprintf( file, ", \\\n '%s.%d.gpt' using %d:%d w l lc rgb '#%02x%02x%02x' title 't%d'", basename, i,
444 [ + - ]: 10 : haxis + 1, vaxis + 1, red( i, count ), green( i, count ), blue( i, count ), i );
445 : : }
446 [ + - ]: 1 : fprintf( file, "\n" );
447 : :
448 [ + - ][ + - ]: 1 : fclose( file );
449 : : }
450 : :
451 : 0 : void write_stl( Mesh* mesh, const char* filename, MsqError& err )
452 : : {
453 : : // Open the file
454 [ # # ][ # # ]: 0 : ofstream file( filename );
455 [ # # ][ # # ]: 0 : if( !file )
456 : : {
457 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
458 : 0 : return;
459 : : }
460 : :
461 : : // Write header
462 : : char header[70];
463 : 0 : sprintf( header, "Mesquite%d", rand() );
464 [ # # ][ # # ]: 0 : file << "solid " << header << endl;
[ # # ]
465 : :
466 [ # # ][ # # ]: 0 : MsqVertex coords[3];
467 [ # # ][ # # ]: 0 : std::vector< Mesh::VertexHandle > verts( 3 );
468 [ # # ][ # # ]: 0 : std::vector< size_t > offsets( 2 );
469 : :
470 : : // Iterate over all elements
471 : 0 : size_t count = 0;
472 [ # # ][ # # ]: 0 : std::vector< Mesh::ElementHandle > elems;
473 : 0 : std::vector< Mesh::ElementHandle >::iterator iter;
474 [ # # ][ # # ]: 0 : mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
475 [ # # ][ # # ]: 0 : for( iter = elems.begin(); iter != elems.end(); ++iter )
[ # # ]
476 : : {
477 : : // Skip non-triangles
478 [ # # ]: 0 : Mesh::ElementHandle elem = *iter;
479 : : EntityTopology type;
480 [ # # ][ # # ]: 0 : mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
481 [ # # ]: 0 : if( type != TRIANGLE ) continue;
482 : 0 : ++count;
483 : :
484 : : // Get vertex coordinates
485 [ # # ][ # # ]: 0 : mesh->elements_get_attached_vertices( &elem, 1, verts, offsets, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
486 [ # # ][ # # ]: 0 : mesh->vertices_get_coordinates( arrptr( verts ), coords, 3, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
487 : :
488 : : // Get triagnle normal
489 [ # # ][ # # ]: 0 : Vector3D n = ( coords[0] - coords[1] ) * ( coords[0] - coords[2] );
[ # # ]
490 [ # # ]: 0 : n.normalize();
491 : :
492 : : // Write triangle
493 [ # # ][ # # ]: 0 : file << "facet normal " << n.x() << " " << n.y() << " " << n.z() << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
494 [ # # ][ # # ]: 0 : file << "outer loop" << endl;
495 [ # # ]: 0 : for( unsigned i = 0; i < 3; ++i )
496 [ # # ][ # # ]: 0 : file << "vertex " << coords[i].x() << " " << coords[i].y() << " " << coords[i].z() << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
497 [ # # ][ # # ]: 0 : file << "endloop" << endl;
498 [ # # ][ # # ]: 0 : file << "endfacet" << endl;
499 : : }
500 : :
501 [ # # ][ # # ]: 0 : file << "endsolid " << header << endl;
[ # # ]
502 : :
503 [ # # ]: 0 : file.close();
504 [ # # ]: 0 : if( count == 0 )
505 : : {
506 : 0 : std::remove( filename );
507 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Mesh contains no triangles", MsqError::INVALID_STATE );
[ # # ]
508 : 0 : }
509 : : }
510 : :
511 : 0 : Projection::Projection( PlanarDomain* domain )
512 : : {
513 [ # # ]: 0 : Vector3D normal;
514 [ # # ]: 0 : domain->vertex_normal_at( 0, normal );
515 [ # # ]: 0 : init( normal );
516 : 0 : }
517 : :
518 : 0 : Projection::Projection( const Vector3D& n )
519 : : {
520 [ # # ]: 0 : init( n );
521 : 0 : }
522 : :
523 : 0 : Projection::Projection( const Vector3D& n, const Vector3D& up )
524 : : {
525 [ # # ]: 0 : init( n, up );
526 : 0 : }
527 : :
528 : 0 : Projection::Projection( Axis h, Axis v )
529 : : {
530 [ # # ][ # # ]: 0 : Vector3D horiz( 0, 0, 0 ), vert( 0, 0, 0 );
531 [ # # ]: 0 : horiz[h] = 1.0;
532 [ # # ]: 0 : vert[v] = 1.0;
533 [ # # ][ # # ]: 0 : init( horiz * vert, vert );
534 : 0 : }
535 : :
536 : 0 : void Projection::init( const Vector3D& n )
537 : : {
538 : : // Choose an "up" direction
539 : :
540 : 0 : Axis max = X;
541 [ # # ]: 0 : for( Axis i = Y; i <= Z; i = ( Axis )( i + 1 ) )
542 [ # # ][ # # ]: 0 : if( fabs( n[i] ) > fabs( n[max] ) ) max = i;
[ # # ]
543 : :
544 : : Axis up;
545 [ # # ]: 0 : if( max == Y )
546 : 0 : up = Z;
547 : : else
548 : 0 : up = Y;
549 : :
550 : : // Initialize rotation matrix
551 : :
552 [ # # ]: 0 : Vector3D up_vect( 0, 0, 0 );
553 [ # # ]: 0 : up_vect[up] = 1.0;
554 [ # # ]: 0 : init( n, up_vect );
555 : 0 : }
556 : :
557 : 0 : void Projection::init( const Vector3D& n1, const Vector3D& up1 )
558 : : {
559 [ # # ]: 0 : MsqError err;
560 [ # # ][ # # ]: 0 : const Vector3D n = n1 / n1.length();
561 [ # # ][ # # ]: 0 : const Vector3D u = up1 / up1.length();
562 : :
563 : : // Rotate for projection
564 [ # # ]: 0 : const Vector3D z( 0., 0., 1. );
565 [ # # ]: 0 : Vector3D r = n * z;
566 [ # # ]: 0 : double angle = r.interior_angle( n, z, err );
567 [ # # ]: 0 : Matrix3D rot1 = rotation( r, angle );
568 : :
569 : : // In-plane rotation for up vector
570 [ # # ][ # # ]: 0 : Vector3D pu = u - n * ( n % u );
[ # # ]
571 [ # # ]: 0 : Vector3D y( 0., 1., 0. );
572 [ # # ]: 0 : angle = z.interior_angle( pu, y, err );
573 [ # # ]: 0 : Matrix3D rot2 = rotation( z, angle );
574 : :
575 [ # # ][ # # ]: 0 : this->myTransform = rot1 * rot2;
576 : 0 : }
577 : :
578 : 0 : Matrix3D Projection::rotation( const Vector3D& axis, double angle )
579 : : {
580 : 0 : const double c = cos( angle );
581 : 0 : const double s = sin( angle );
582 [ # # ]: 0 : const double x = axis.x();
583 [ # # ]: 0 : const double y = axis.y();
584 [ # # ]: 0 : const double z = axis.z();
585 : :
586 : 0 : const double xform[9] = { c + x * x * ( 1 - c ), -z * s + x * y * ( 1 - c ), y * s + x * z * ( 1 - c ),
587 : 0 : z * s + x * y * ( 1 - c ), c + y * y * ( 1 - c ), -x * s + y * z * ( 1 - c ),
588 : 0 : -y * s + x * z * ( 1 - c ), x * s + y * z * ( 1 - c ), c + z * z * ( 1 - c ) };
589 [ # # ]: 0 : return Matrix3D( xform );
590 : : }
591 : :
592 : 0 : void Projection::project( const Vector3D& p, float& h, float& v )
593 : : {
594 [ # # ]: 0 : Vector3D pr = myTransform * p;
595 [ # # ]: 0 : h = (float)pr.x();
596 [ # # ]: 0 : v = (float)pr.y();
597 : 0 : }
598 : :
599 : 0 : Transform2D::Transform2D( PatchData* pd, Projection& projection, unsigned width, unsigned height, bool flip )
600 [ # # ]: 0 : : myProj( projection ), vertSign( flip ? -1 : 1 )
601 : : {
602 : : // Get the bounding box of the projected points
603 : : float w_max, w_min, h_max, h_min;
604 : 0 : w_max = h_max = -std::numeric_limits< float >::max();
605 : 0 : w_min = h_min = std::numeric_limits< float >::max();
606 [ # # ]: 0 : MsqError err;
607 [ # # ]: 0 : const MsqVertex* verts = pd->get_vertex_array( err );
608 [ # # ]: 0 : const size_t num_vert = pd->num_nodes();
609 [ # # ]: 0 : for( unsigned i = 0; i < num_vert; ++i )
610 : : {
611 : : float w, h;
612 [ # # ]: 0 : myProj.project( verts[i], w, h );
613 [ # # ]: 0 : if( w > w_max ) w_max = w;
614 [ # # ]: 0 : if( w < w_min ) w_min = w;
615 [ # # ]: 0 : if( h > h_max ) h_max = h;
616 [ # # ]: 0 : if( h < h_min ) h_min = h;
617 : : }
618 : :
619 : : // Determine the scale factor
620 : 0 : const float w_scale = (float)width / ( w_max - w_min );
621 : 0 : const float h_scale = (float)height / ( h_max - h_min );
622 [ # # ]: 0 : myScale = w_scale > h_scale ? h_scale : w_scale;
623 : :
624 : : // Determine offset
625 : 0 : horizOffset = -(int)( myScale * w_min );
626 [ # # ]: 0 : vertOffset = -(int)( myScale * ( flip ? -h_max : h_min ) );
627 : :
628 : : // Determine bounding box
629 : 0 : horizMax = (int)( w_max * myScale ) + horizOffset;
630 [ # # ]: 0 : vertMax = (int)( ( flip ? -h_min : h_max ) * myScale ) + vertOffset;
631 : 0 : }
632 : :
633 : 0 : Transform2D::Transform2D( const Vector3D* verts, size_t num_vert, Projection& projection, unsigned width,
634 : : unsigned height )
635 : 0 : : myProj( projection ), vertSign( 1 )
636 : : {
637 : : // Get the bounding box of the projected points
638 : : float w_max, w_min, h_max, h_min;
639 : 0 : w_max = h_max = -std::numeric_limits< float >::max();
640 : 0 : w_min = h_min = std::numeric_limits< float >::max();
641 [ # # ]: 0 : for( unsigned i = 0; i < num_vert; ++i )
642 : : {
643 : : float w, h;
644 [ # # ]: 0 : myProj.project( verts[i], w, h );
645 [ # # ]: 0 : if( w > w_max ) w_max = w;
646 [ # # ]: 0 : if( w < w_min ) w_min = w;
647 [ # # ]: 0 : if( h > h_max ) h_max = h;
648 [ # # ]: 0 : if( h < h_min ) h_min = h;
649 : : }
650 : :
651 : : // Determine the scale factor
652 : 0 : const float w_scale = (float)width / ( w_max - w_min );
653 : 0 : const float h_scale = (float)height / ( h_max - h_min );
654 [ # # ]: 0 : myScale = w_scale > h_scale ? h_scale : w_scale;
655 : :
656 : : // Determine offset
657 : 0 : horizOffset = -(int)( myScale * w_min );
658 : 0 : vertOffset = -(int)( myScale * h_min );
659 : :
660 : : // Determine bounding box
661 : 0 : horizMax = (int)( w_max * myScale ) + horizOffset;
662 : 0 : vertMax = (int)( h_max * myScale ) + vertOffset;
663 : 0 : }
664 : :
665 : 0 : void Transform2D::transform( const Vector3D& coords, int& horizontal, int& vertical ) const
666 : : {
667 : : float horiz, vert;
668 [ # # ]: 0 : myProj.project( coords, horiz, vert );
669 : 0 : horizontal = (int)( myScale * horiz ) + horizOffset;
670 : 0 : vertical = vertSign * (int)( myScale * vert ) + vertOffset;
671 : 0 : }
672 : :
673 : : /** Write quadratic edge shape in PostScript format.
674 : : *
675 : : * Given the three points composing a quadratic mesh edge,
676 : : * write the cubic Bezier curve of the same shape in
677 : : * PostScript format. The formulas for P1 and P2
678 : : * at the start of this function will result in the cubic
679 : : * terms of the Bezier curve dropping out, leaving the
680 : : * quadratic curve matching the edge shape function as
681 : : * described in Section 3.6 of Hughes. (If you're attempting
682 : : * to verify this, don't forget to adjust for the different
683 : : * parameter ranges: \f$ \xi = 2 t - 1 \f$).
684 : : */
685 : 0 : static void write_eps_quadratic_edge( ostream& s, Transform2D& xform, Vector3D start, Vector3D mid, Vector3D end )
686 : : {
687 [ # # ][ # # ]: 0 : Vector3D P1 = 1. / 3 * ( 4 * mid - end );
[ # # ]
688 [ # # ][ # # ]: 0 : Vector3D P2 = 1. / 3 * ( 4 * mid - start );
[ # # ]
689 : :
690 : : int x, y;
691 [ # # ]: 0 : xform.transform( start, x, y );
692 [ # # ][ # # ]: 0 : s << x << ' ' << y << " moveto" << endl;
[ # # ][ # # ]
[ # # ]
693 [ # # ]: 0 : xform.transform( P1, x, y );
694 [ # # ][ # # ]: 0 : s << x << ' ' << y << ' ';
[ # # ][ # # ]
695 [ # # ]: 0 : xform.transform( P2, x, y );
696 [ # # ][ # # ]: 0 : s << x << ' ' << y << ' ';
[ # # ][ # # ]
697 [ # # ]: 0 : xform.transform( end, x, y );
698 [ # # ][ # # ]: 0 : s << x << ' ' << y << " curveto" << endl;
[ # # ][ # # ]
[ # # ]
699 : 0 : }
700 : :
701 : 0 : void write_eps( Mesh* mesh, const char* filename, Projection proj, MsqError& err, int width, int height )
702 : : {
703 : : // Get a global patch
704 [ # # ]: 0 : PatchData pd;
705 [ # # ]: 0 : pd.set_mesh( mesh );
706 [ # # ][ # # ]: 0 : pd.fill_global_patch( err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
707 : :
708 [ # # ]: 0 : Transform2D transf( &pd, proj, width, height, false );
709 : :
710 : : // Open the file
711 [ # # ][ # # ]: 0 : ofstream s( filename );
[ # # ]
712 [ # # ][ # # ]: 0 : if( !s )
713 : : {
714 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
715 : 0 : return;
716 : : }
717 : :
718 : : // Write header
719 [ # # ][ # # ]: 0 : s << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
720 [ # # ][ # # ]: 0 : s << "%%Creator: Mesquite" << endl;
721 [ # # ][ # # ]: 0 : s << "%%Title: Mesquite " << endl;
722 [ # # ][ # # ]: 0 : s << "%%DocumentData: Clean7Bit" << endl;
723 [ # # ][ # # ]: 0 : s << "%%Origin: 0 0" << endl;
724 [ # # ][ # # ]: 0 : s << "%%BoundingBox: 0 0 " << transf.max_horizontal() << ' ' << transf.max_vertical() << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
725 [ # # ][ # # ]: 0 : s << "%%Pages: 1" << endl;
726 : :
727 [ # # ][ # # ]: 0 : s << "%%BeginProlog" << endl;
728 [ # # ][ # # ]: 0 : s << "save" << endl;
729 [ # # ][ # # ]: 0 : s << "countdictstack" << endl;
730 [ # # ][ # # ]: 0 : s << "mark" << endl;
731 [ # # ][ # # ]: 0 : s << "newpath" << endl;
732 [ # # ][ # # ]: 0 : s << "/showpage {} def" << endl;
733 [ # # ][ # # ]: 0 : s << "/setpagedevice {pop} def" << endl;
734 [ # # ][ # # ]: 0 : s << "%%EndProlog" << endl;
735 : :
736 [ # # ][ # # ]: 0 : s << "%%Page: 1 1" << endl;
737 [ # # ][ # # ]: 0 : s << "1 setlinewidth" << endl;
738 [ # # ][ # # ]: 0 : s << "0.0 setgray" << endl;
739 : :
740 : : // Write mesh edges
741 [ # # ][ # # ]: 0 : EdgeIterator iter( &pd, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
742 [ # # ][ # # ]: 0 : while( !iter.is_at_end() )
743 : : {
744 : : int s_w, s_h, e_w, e_h;
745 [ # # ][ # # ]: 0 : transf.transform( iter.start(), s_w, s_h );
746 [ # # ][ # # ]: 0 : transf.transform( iter.end(), e_w, e_h );
747 : :
748 [ # # ][ # # ]: 0 : s << "newpath" << endl;
749 [ # # ][ # # ]: 0 : s << s_w << ' ' << s_h << " moveto" << endl;
[ # # ][ # # ]
[ # # ]
750 : :
751 [ # # ][ # # ]: 0 : if( !iter.mid() ) { s << e_w << ' ' << e_h << " lineto" << endl; }
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
752 : : else
753 : : {
754 [ # # ][ # # ]: 0 : write_eps_quadratic_edge( s, transf, iter.start(), *iter.mid(), iter.end() );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
755 : : // draw rings at mid-edge node location
756 : : // transf.transform( *(iter.mid()), w1, h1 );
757 : : // s << w1+2 << ' ' << h1 << " moveto" << endl;
758 : : // s << w1 << ' ' << h1 << " 2 0 360 arc" << endl;
759 : : }
760 [ # # ][ # # ]: 0 : s << "stroke" << endl;
761 : :
762 [ # # ][ # # ]: 0 : iter.step( err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
763 : : }
764 : :
765 : : // Write footer
766 [ # # ][ # # ]: 0 : s << "%%Trailer" << endl;
767 [ # # ][ # # ]: 0 : s << "cleartomark" << endl;
768 [ # # ][ # # ]: 0 : s << "countdictstack" << endl;
769 [ # # ][ # # ]: 0 : s << "exch sub { end } repeat" << endl;
770 [ # # ][ # # ]: 0 : s << "restore" << endl;
771 [ # # ][ # # ]: 0 : s << "%%EOF" << endl;
[ # # ]
772 : : }
773 : :
774 : : /** Quadratic triangle shape function for use in write_eps_triangle */
775 : 0 : static double tN0( double r, double s )
776 : : {
777 : 0 : double t = 1 - r - s;
778 : 0 : return t * ( 2 * t - 1 );
779 : : }
780 : 0 : static double tN1( double r, double )
781 : : {
782 : 0 : return r * ( 2 * r - 1 );
783 : : }
784 : 0 : static double tN2( double, double s )
785 : : {
786 : 0 : return s * ( 2 * s - 1 );
787 : : }
788 : 0 : static double tN3( double r, double s )
789 : : {
790 : 0 : double t = 1 - r - s;
791 : 0 : return 4 * r * t;
792 : : }
793 : 0 : static double tN4( double r, double s )
794 : : {
795 : 0 : return 4 * r * s;
796 : : }
797 : 0 : static double tN5( double r, double s )
798 : : {
799 : 0 : double t = 1 - r - s;
800 : 0 : return 4 * s * t;
801 : : }
802 : : /** Quadratic triangle shape function for use in write_eps_triangle */
803 : 0 : static Vector3D quad_tri_pt( double r, double s, const Vector3D* coords )
804 : : {
805 : 0 : Vector3D result = tN0( r, s ) * coords[0];
806 [ # # ]: 0 : result += tN1( r, s ) * coords[1];
807 [ # # ]: 0 : result += tN2( r, s ) * coords[2];
808 [ # # ]: 0 : result += tN3( r, s ) * coords[3];
809 [ # # ]: 0 : result += tN4( r, s ) * coords[4];
810 [ # # ]: 0 : result += tN5( r, s ) * coords[5];
811 : 0 : return result;
812 : : }
813 : :
814 : 0 : void write_eps_triangle( Mesh* mesh, Mesh::ElementHandle elem, const char* filename, bool draw_iso_lines,
815 : : bool draw_nodes, MsqError& err, int width, int height )
816 : : {
817 : : // Get triangle vertices
818 [ # # ][ # # ]: 0 : MsqVertex coords[6];
819 : : EntityTopology type;
820 [ # # ][ # # ]: 0 : mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
821 [ # # ]: 0 : if( type != TRIANGLE )
822 : : {
823 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Invalid element type", MsqError::UNSUPPORTED_ELEMENT );
824 : 0 : return;
825 : : }
826 [ # # ]: 0 : std::vector< Mesh::VertexHandle > verts;
827 [ # # ][ # # ]: 0 : std::vector< size_t > junk;
828 [ # # ][ # # ]: 0 : mesh->elements_get_attached_vertices( &elem, 1, verts, junk, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
829 [ # # ][ # # ]: 0 : if( verts.size() != 3 && verts.size() != 6 )
[ # # ]
830 : : {
831 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( "Invalid element type", MsqError::UNSUPPORTED_ELEMENT );
832 : 0 : return;
833 : : }
834 [ # # ][ # # ]: 0 : mesh->vertices_get_coordinates( arrptr( verts ), coords, verts.size(), err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
835 : :
836 [ # # ][ # # ]: 0 : Vector3D coords2[6];
837 [ # # ]: 0 : std::copy( coords, coords + verts.size(), coords2 );
838 : :
839 [ # # ][ # # ]: 0 : std::vector< bool > fixed( verts.size(), false );
840 [ # # ]: 0 : if( draw_nodes )
841 : : {
842 [ # # ][ # # ]: 0 : mesh->vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
843 : : }
844 [ # # ][ # # ]: 0 : write_eps_triangle( coords2, verts.size(), filename, draw_iso_lines, draw_nodes, err, fixed, width, height );
845 : : }
846 : :
847 : 0 : void write_eps_triangle( const Vector3D* coords, size_t num_vtx, const char* filename, bool draw_iso_lines,
848 : : bool draw_nodes, MsqError& err, const std::vector< bool >& fixed, int width, int height )
849 : : {
850 : 0 : const int PT_RAD = 3; // radius of circles for drawing nodes, in points
851 : 0 : const int PAD = PT_RAD + 2; // margin in points
852 : 0 : const double EDGE_GRAY = 0.0; // color for triangle edges, 0.0 => black
853 : 0 : const double ISO_GRAY = 0.7; // color for parameter iso-lines
854 : 0 : const double NODE_GRAY = 0.0; // color for node circle
855 : 0 : const double FIXED_GRAY = 1.0; // color to fill fixed nodes with, 1.0 => white
856 : 0 : const double FREE_GRAY = 0.0; // color to fill free nodes with
857 : :
858 [ # # ]: 0 : Projection proj( X, Y );
859 [ # # ]: 0 : Transform2D transf( coords, num_vtx, proj, width, height );
860 : :
861 : : // Open the file
862 [ # # ][ # # ]: 0 : ofstream str( filename );
[ # # ]
863 [ # # ][ # # ]: 0 : if( !str )
864 : : {
865 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
866 : 0 : return;
867 : : }
868 : :
869 : : // Write header
870 [ # # ][ # # ]: 0 : str << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
871 [ # # ][ # # ]: 0 : str << "%%Creator: Mesquite" << endl;
872 [ # # ][ # # ]: 0 : str << "%%Title: Mesquite " << endl;
873 [ # # ][ # # ]: 0 : str << "%%DocumentData: Clean7Bit" << endl;
874 [ # # ][ # # ]: 0 : str << "%%Origin: 0 0" << endl;
875 [ # # ][ # # ]: 0 : str << "%%BoundingBox: " << -PAD << ' ' << -PAD << ' ' << transf.max_horizontal() + PAD << ' '
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ][ # # ]
876 [ # # ][ # # ]: 0 : << transf.max_vertical() + PAD << endl;
[ # # ]
877 [ # # ][ # # ]: 0 : str << "%%Pages: 1" << endl;
878 : :
879 [ # # ][ # # ]: 0 : str << "%%BeginProlog" << endl;
880 [ # # ][ # # ]: 0 : str << "save" << endl;
881 [ # # ][ # # ]: 0 : str << "countdictstack" << endl;
882 [ # # ][ # # ]: 0 : str << "mark" << endl;
883 [ # # ][ # # ]: 0 : str << "newpath" << endl;
884 [ # # ][ # # ]: 0 : str << "/showpage {} def" << endl;
885 [ # # ][ # # ]: 0 : str << "/setpagedevice {pop} def" << endl;
886 [ # # ][ # # ]: 0 : str << "%%EndProlog" << endl;
887 : :
888 [ # # ][ # # ]: 0 : str << "%%Page: 1 1" << endl;
889 [ # # ][ # # ]: 0 : str << "1 setlinewidth" << endl;
890 [ # # ][ # # ]: 0 : str << EDGE_GRAY << " setgray" << endl;
[ # # ]
891 : :
892 : 0 : const double h = 0.5, t = 1.0 / 3.0, w = 2.0 / 3.0, s = 1. / 6, f = 5. / 6;
893 : 0 : const int NUM_ISO = 15;
894 : : const double iso_params[NUM_ISO][2][2] = {
895 : : { { h, 0 }, { h, h } }, // r = 1/2
896 : : { { t, 0 }, { t, w } }, // r = 1/3
897 : : { { w, 0 }, { w, t } }, // r = 2/3
898 : : { { s, 0 }, { s, f } }, // r = 1/6
899 : : { { f, 0 }, { f, s } }, // r = 5/6
900 : : { { 0, h }, { h, h } }, // s = 1/2
901 : : { { 0, t }, { w, t } }, // s = 1/3
902 : : { { 0, w }, { t, w } }, // s = 2/3
903 : : { { 0, s }, { f, s } }, // s = 1/6
904 : : { { 0, f }, { s, f } }, // s = 5/6
905 : : { { 0, h }, { h, 0 } }, // t = 1 - r - s = 1/2
906 : : { { 0, w }, { w, 0 } }, // t = 1 - r - s = 1/3
907 : : { { 0, t }, { t, 0 } }, // t = 1 - r - s = 2/3
908 : : { { 0, f }, { f, 0 } }, // t = 1 - r - s = 1/6
909 : : { { 0, s }, { s, 0 } } // t = 1 - r - s = 5/6
910 : 0 : };
911 : :
912 [ # # ]: 0 : if( num_vtx == 3 )
913 : : {
914 : : int x[3], y[3];
915 [ # # ]: 0 : for( size_t i = 0; i < 3; ++i )
916 [ # # ]: 0 : transf.transform( coords[i], x[i], y[i] );
917 : :
918 [ # # ][ # # ]: 0 : str << "newpath" << endl;
919 [ # # ][ # # ]: 0 : str << x[0] << ' ' << y[0] << " moveto" << endl;
[ # # ][ # # ]
[ # # ]
920 [ # # ][ # # ]: 0 : str << x[1] << ' ' << y[1] << " lineto" << endl;
[ # # ][ # # ]
[ # # ]
921 [ # # ][ # # ]: 0 : str << x[2] << ' ' << y[2] << " lineto" << endl;
[ # # ][ # # ]
[ # # ]
922 [ # # ][ # # ]: 0 : str << x[0] << ' ' << y[0] << " lineto" << endl;
[ # # ][ # # ]
[ # # ]
923 [ # # ][ # # ]: 0 : str << "stroke" << endl;
924 : :
925 [ # # ]: 0 : if( draw_iso_lines )
926 : : {
927 [ # # ][ # # ]: 0 : str << ISO_GRAY << " setgray" << endl;
[ # # ]
928 [ # # ][ # # ]: 0 : str << "newpath" << endl;
929 [ # # ]: 0 : for( int i = 0; i < NUM_ISO; ++i )
930 : : {
931 : 0 : double R[2] = { iso_params[i][0][0], iso_params[i][1][0] };
932 : 0 : double S[2] = { iso_params[i][0][1], iso_params[i][1][1] };
933 : 0 : double T[2] = { 1 - R[0] - S[0], 1 - R[1] - S[1] };
934 [ # # ][ # # ]: 0 : Vector3D p[2] = { T[0] * coords[0] + R[0] * coords[1] + S[0] * coords[2],
[ # # ][ # # ]
935 [ # # ][ # # ]: 0 : T[1] * coords[0] + R[1] * coords[1] + S[1] * coords[2] };
[ # # ][ # # ]
[ # # ][ # # ]
936 [ # # ]: 0 : transf.transform( p[0], x[0], y[0] );
937 [ # # ]: 0 : transf.transform( p[1], x[1], y[1] );
938 [ # # ][ # # ]: 0 : str << x[0] << ' ' << y[0] << " moveto" << endl;
[ # # ][ # # ]
[ # # ]
939 [ # # ][ # # ]: 0 : str << x[1] << ' ' << y[1] << " lineto" << endl;
[ # # ][ # # ]
[ # # ]
940 : : }
941 [ # # ][ # # ]: 0 : str << " stroke" << endl;
942 : : }
943 : : }
944 [ # # ]: 0 : else if( num_vtx == 6 )
945 : : {
946 [ # # ][ # # ]: 0 : str << "newpath" << endl;
947 [ # # ][ # # ]: 0 : write_eps_quadratic_edge( str, transf, coords[0], coords[3], coords[1] );
[ # # ][ # # ]
948 [ # # ][ # # ]: 0 : write_eps_quadratic_edge( str, transf, coords[1], coords[4], coords[2] );
[ # # ][ # # ]
949 [ # # ][ # # ]: 0 : write_eps_quadratic_edge( str, transf, coords[2], coords[5], coords[0] );
[ # # ][ # # ]
950 [ # # ][ # # ]: 0 : str << "stroke" << endl;
951 : :
952 [ # # ]: 0 : if( draw_iso_lines )
953 : : {
954 [ # # ][ # # ]: 0 : str << ISO_GRAY << " setgray" << endl;
[ # # ]
955 [ # # ][ # # ]: 0 : str << "newpath" << endl;
956 [ # # ]: 0 : for( int i = 0; i < NUM_ISO; ++i )
957 : : {
958 : 0 : double R[3] = { iso_params[i][0][0], 0, iso_params[i][1][0] };
959 : 0 : double S[3] = { iso_params[i][0][1], 0, iso_params[i][1][1] };
960 : 0 : R[1] = 0.5 * ( R[0] + R[2] );
961 : 0 : S[1] = 0.5 * ( S[0] + S[2] );
962 : : Vector3D p[3] = { quad_tri_pt( R[0], S[0], coords ), quad_tri_pt( R[1], S[1], coords ),
963 [ # # ][ # # ]: 0 : quad_tri_pt( R[2], S[2], coords ) };
[ # # ]
964 [ # # ][ # # ]: 0 : write_eps_quadratic_edge( str, transf, p[0], p[1], p[2] );
[ # # ][ # # ]
965 : : }
966 [ # # ][ # # ]: 0 : str << " stroke" << endl;
967 : : }
968 : : }
969 : :
970 [ # # ]: 0 : if( draw_nodes )
971 : : {
972 [ # # ]: 0 : for( size_t i = 0; i < num_vtx; ++i )
973 : : {
974 : : int iw, ih;
975 : : // fill interior with either white or black depending
976 : : // on whether or not the vertex is fixed.
977 [ # # ][ # # ]: 0 : if( fixed[i] )
978 [ # # ][ # # ]: 0 : str << FIXED_GRAY << " setgray" << endl;
[ # # ]
979 : : else
980 [ # # ][ # # ]: 0 : str << FREE_GRAY << " setgray" << endl;
[ # # ]
981 [ # # ]: 0 : transf.transform( coords[i], iw, ih );
982 [ # # ][ # # ]: 0 : str << w + PT_RAD << ' ' << ih << " moveto" << endl;
[ # # ][ # # ]
[ # # ]
983 [ # # ][ # # ]: 0 : str << iw << ' ' << ih << ' ' << PT_RAD << " 0 360 arc" << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
984 [ # # ][ # # ]: 0 : str << "closepath fill" << endl;
985 [ # # ][ # # ]: 0 : str << NODE_GRAY << " setgray" << endl;
[ # # ]
986 [ # # ][ # # ]: 0 : str << "newpath" << endl;
987 [ # # ][ # # ]: 0 : str << iw << ' ' << ih << ' ' << PT_RAD << " 0 360 arc" << endl;
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
988 [ # # ][ # # ]: 0 : str << "stroke" << endl;
989 : : }
990 : : }
991 : :
992 : : // Write footer
993 [ # # ][ # # ]: 0 : str << "%%Trailer" << endl;
994 [ # # ][ # # ]: 0 : str << "cleartomark" << endl;
995 [ # # ][ # # ]: 0 : str << "countdictstack" << endl;
996 [ # # ][ # # ]: 0 : str << "exch sub { end } repeat" << endl;
997 [ # # ][ # # ]: 0 : str << "restore" << endl;
998 [ # # ][ # # ]: 0 : str << "%%EOF" << endl;
[ # # ]
999 : : }
1000 : :
1001 : 0 : void write_svg( Mesh* mesh, const char* filename, Projection proj, MsqError& err )
1002 : : {
1003 : : // Get a global patch
1004 [ # # ]: 0 : PatchData pd;
1005 [ # # ]: 0 : pd.set_mesh( mesh );
1006 [ # # ][ # # ]: 0 : pd.fill_global_patch( err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
1007 : :
1008 [ # # ]: 0 : Transform2D transf( &pd, proj, 400, 400, true );
1009 : :
1010 : : // Open the file
1011 [ # # ][ # # ]: 0 : ofstream file( filename );
[ # # ]
1012 [ # # ][ # # ]: 0 : if( !file )
1013 : : {
1014 [ # # ][ # # ]: 0 : MSQ_SETERR( err )( MsqError::FILE_ACCESS );
1015 : 0 : return;
1016 : : }
1017 : :
1018 : : // Write header
1019 [ # # ][ # # ]: 0 : file << "<?xml version=\"1.0\" standalone=\"no\"?>" << endl;
1020 [ # # ]: 0 : file << "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" "
1021 [ # # ][ # # ]: 0 : << "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << endl;
1022 [ # # ]: 0 : file << endl;
1023 [ # # ][ # # ]: 0 : file << "<svg width=\"" << transf.max_horizontal() << "\" height=\"" << transf.max_vertical()
[ # # ][ # # ]
[ # # ][ # # ]
1024 [ # # ][ # # ]: 0 : << "\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" << endl;
1025 : :
1026 : : // Write mesh edges
1027 [ # # ][ # # ]: 0 : EdgeIterator iter( &pd, err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
[ # # ]
1028 [ # # ][ # # ]: 0 : while( !iter.is_at_end() )
1029 : : {
1030 : : int s_w, s_h, e_w, e_h;
1031 [ # # ][ # # ]: 0 : transf.transform( iter.start(), s_w, s_h );
1032 [ # # ][ # # ]: 0 : transf.transform( iter.end(), e_w, e_h );
1033 : :
1034 [ # # ]: 0 : file << "<line "
1035 [ # # ][ # # ]: 0 : << "x1=\"" << s_w << "\" "
[ # # ]
1036 [ # # ][ # # ]: 0 : << "y1=\"" << s_h << "\" "
[ # # ]
1037 [ # # ][ # # ]: 0 : << "x2=\"" << e_w << "\" "
[ # # ]
1038 [ # # ][ # # ]: 0 : << "y2=\"" << e_h << "\" "
[ # # ]
1039 [ # # ]: 0 : << " style=\"stroke:rgb(99,99,99);stroke-width:2\""
1040 [ # # ][ # # ]: 0 : << "/>" << endl;
1041 : :
1042 [ # # ][ # # ]: 0 : iter.step( err );MSQ_ERRRTN( err );
[ # # ][ # # ]
[ # # ][ # # ]
1043 : : }
1044 : :
1045 : : // Write footer
1046 [ # # ][ # # ]: 0 : file << "</svg>" << endl;
[ # # ]
1047 : : }
1048 : :
1049 : : } // namespace MeshWriter
1050 : :
1051 [ + - ][ + - ]: 120 : } // namespace MBMesquite
1052 : :
1053 : : #endif
|
