MOAB: Mesh Oriented datABase  (version 5.4.0)
MeshWriter.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2005 Lawrence Livermore National Laboratory.  Under
00005     the terms of Contract B545069 with the University of Wisconsin --
00006     Madison, Lawrence Livermore National Laboratory retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 #ifndef MSQ_MESH_WRITER_CPP
00028 #define MSQ_MESH_WRITER_CPP
00029 
00030 #include "MeshWriter.hpp"
00031 #include "Mesquite.hpp"
00032 #include "MeshImpl.hpp"
00033 #include "MsqError.hpp"
00034 #include "PatchData.hpp"
00035 #include "PlanarDomain.hpp"
00036 #include "VtkTypeInfo.hpp"
00037 #include "EdgeIterator.hpp"
00038 
00039 #include <memory>
00040 #include <limits>
00041 #include <vector>
00042 #include <algorithm>
00043 #include <fstream>
00044 #include <string>
00045 #include <iomanip>
00046 using namespace std;
00047 
00048 #include <cstdio>
00049 
00050 namespace MBMesquite
00051 {
00052 
00053 namespace MeshWriter
00054 {
00055 
00056     /**\brief Transform from coordinates in the XY-plane
00057      *        to graphics coordinates.
00058      */
00059     class Transform2D
00060     {
00061       public:
00062         Transform2D( PatchData* pd, Projection& proj, unsigned width, unsigned height, bool flip_about_horizontal );
00063 
00064         Transform2D( const Vector3D* verts, size_t num_vert, Projection& projection, unsigned width, unsigned height );
00065 
00066         void transform( const Vector3D& coords, int& horizontal, int& vertical ) const;
00067 
00068         int max_horizontal() const
00069         {
00070             return horizMax;
00071         }
00072         int max_vertical() const
00073         {
00074             return vertMax;
00075         }
00076 
00077       private:
00078         Projection& myProj;
00079         float myScale;
00080         int horizOffset, vertOffset;
00081         int horizMax, vertMax;
00082         int vertSign;
00083     };
00084 
00085     /* Write VTK file
00086      *
00087      * Copied from src/Mesh/MeshSet.cpp and adapted for removal of
00088      * MeshSet class by J.Kraftcheck on 2005-7-28.
00089      *
00090      * This code is provided mainly for debugging.  A more efficient
00091      * and complete writer implementation is provided in the MeshImpl
00092      * class for saving meshes that were read from a file initially.
00093      */
00094     void write_vtk( Mesh* mesh, const char* out_filename, MsqError& err )
00095     {
00096         if( MeshImpl* msq_mesh = dynamic_cast< MeshImpl* >( mesh ) )
00097         {
00098             msq_mesh->write_vtk( out_filename, err );MSQ_CHKERR( err );
00099             return;
00100         }
00101 
00102         // loads a global patch
00103         PatchData pd;
00104         pd.set_mesh( mesh );
00105         pd.fill_global_patch( err );MSQ_ERRRTN( err );
00106 
00107         // write mesh
00108         write_vtk( pd, out_filename, err );MSQ_CHKERR( err );
00109     }
00110 
00111     void write_vtk( PatchData& pd, const char* out_filename, MsqError& err, const Vector3D* OF_gradient )
00112     {
00113         // Open the file
00114         std::ofstream file( out_filename );
00115         if( !file )
00116         {
00117             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00118             return;
00119         }
00120 
00121         // Write a header
00122         file << "# vtk DataFile Version 2.0\n";
00123         file << "Mesquite Mesh " << out_filename << " .\n";
00124         file << "ASCII\n";
00125         file << "DATASET UNSTRUCTURED_GRID\n";
00126 
00127         // Write vertex coordinates
00128         file << "POINTS " << pd.num_nodes() << " float\n";
00129         size_t i;
00130         for( i = 0; i < pd.num_nodes(); i++ )
00131         {
00132             file << pd.vertex_by_index( i )[0] << ' ' << pd.vertex_by_index( i )[1] << ' ' << pd.vertex_by_index( i )[2]
00133                  << '\n';
00134         }
00135 
00136         // Write out the connectivity table
00137         size_t connectivity_size = 0;
00138         for( i = 0; i < pd.num_elements(); ++i )
00139             connectivity_size += pd.element_by_index( i ).node_count() + 1;
00140 
00141         file << "CELLS " << pd.num_elements() << ' ' << connectivity_size << '\n';
00142         for( i = 0; i < pd.num_elements(); i++ )
00143         {
00144             std::vector< size_t > vtx_indices;
00145             pd.element_by_index( i ).get_node_indices( vtx_indices );
00146 
00147             // Convert native to VTK node order, if not the same
00148             const VtkTypeInfo* info =
00149                 VtkTypeInfo::find_type( pd.element_by_index( i ).get_element_type(), vtx_indices.size(), err );MSQ_ERRRTN( err );
00150             info->mesquiteToVtkOrder( vtx_indices );
00151 
00152             file << vtx_indices.size();
00153             for( std::size_t j = 0; j < vtx_indices.size(); ++j )
00154             {
00155                 file << ' ' << vtx_indices[j];
00156             }
00157             file << '\n';
00158         }
00159 
00160         // Write out the element types
00161         file << "CELL_TYPES " << pd.num_elements() << '\n';
00162         for( i = 0; i < pd.num_elements(); i++ )
00163         {
00164             const VtkTypeInfo* info = VtkTypeInfo::find_type( pd.element_by_index( i ).get_element_type(),
00165                                                               pd.element_by_index( i ).node_count(), err );MSQ_ERRRTN( err );
00166             file << info->vtkType << '\n';
00167         }
00168 
00169         // Write out which points are fixed.
00170         file << "POINT_DATA " << pd.num_nodes() << "\nSCALARS fixed int\nLOOKUP_TABLE default\n";
00171         for( i = 0; i < pd.num_nodes(); ++i )
00172         {
00173             if( pd.vertex_by_index( i ).get_flags() & MsqVertex::MSQ_CULLED )
00174                 file << "1\n";
00175             else
00176                 file << "0\n";
00177         }
00178         file << "SCALARS culled short\nLOOKUP_TABLE default\n";
00179         for( i = 0; i < pd.num_nodes(); ++i )
00180         {
00181             if( pd.vertex_by_index( i ).is_free_vertex() )
00182                 file << "0\n";
00183             else
00184                 file << "1\n";
00185         }
00186 
00187         if( OF_gradient )
00188         {
00189             file << "VECTORS gradient double\n";
00190             for( i = 0; i < pd.num_free_vertices(); ++i )
00191                 file << OF_gradient[i].x() << " " << OF_gradient[i].y() << " " << OF_gradient[i].z() << "\n";
00192             for( i = pd.num_free_vertices(); i < pd.num_nodes(); ++i )
00193                 file << "0.0 0.0 0.0\n";
00194         }
00195 
00196         // Close the file
00197         file.close();
00198     }
00199 
00200     void write_gnuplot( Mesh* mesh, const char* out_filebase, MsqError& err )
00201     {
00202         // loads a global patch
00203         PatchData pd;
00204         pd.set_mesh( mesh );
00205         pd.fill_global_patch( err );MSQ_ERRRTN( err );
00206         write_gnuplot( pd, out_filebase, err );
00207     }
00208 
00209     void write_gnuplot( Mesh* mesh, std::vector< Mesh::ElementHandle >& elems, const char* out_filebase, MsqError& err )
00210     {
00211         // loads a global patch
00212         PatchData pd;
00213         pd.set_mesh( mesh );
00214         std::vector< Mesh::VertexHandle > verts;
00215         pd.set_mesh_entities( elems, verts, err );MSQ_ERRRTN( err );
00216         write_gnuplot( pd, out_filebase, err );
00217     }
00218 
00219     /*  Writes a gnuplot file directly from the MeshSet.
00220      *  This means that any mesh imported successfully into Mesquite
00221      *  can be outputed in gnuplot format.
00222      *
00223      *  Within gnuplot, use \b plot 'file1.gpt' w l, 'file2.gpt' w l
00224      *
00225      *  This is not geared for performance, since it has to load a global Patch from
00226      *  the mesh to write a mesh file.
00227      *
00228      * Copied from src/Mesh/MeshSet.cpp and adapted for removal of
00229      * MeshSet class by J.Kraftcheck on 2005-7-28.
00230      *
00231      * Re-written to use EdgeIterator by J.Kraftcheck on 2009-6-11
00232      */
00233     void write_gnuplot( PatchData& pd, const char* out_filebase, MsqError& err )
00234     {
00235         // Open the file
00236         std::string out_filename = out_filebase;
00237         out_filename += ".gpt";
00238         std::ofstream file( out_filename.c_str() );
00239         if( !file )
00240         {
00241             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00242             return;
00243         }
00244 
00245         EdgeIterator edges( &pd, err );MSQ_ERRRTN( err );
00246 
00247         // Write a header
00248         file << "\n";
00249 
00250         while( !edges.is_at_end() )
00251         {
00252             const Vector3D& s = edges.start();
00253             const Vector3D& e = edges.end();
00254             const Vector3D* m = edges.mid();
00255 
00256             file << s[0] << ' ' << s[1] << ' ' << s[2] << std::endl;
00257             if( m ) file << ( *m )[0] << ' ' << ( *m )[1] << ' ' << ( *m )[2] << std::endl;
00258             file << e[0] << ' ' << e[1] << ' ' << e[2] << std::endl;
00259             file << std::endl << std::endl;
00260 
00261             edges.step( err );MSQ_ERRRTN( err );
00262         }
00263 
00264         // Close the file
00265         file.close();
00266     }
00267 
00268     /** Helper function for write_gnuplot_animator and write_gnuplot_overlay
00269      *
00270      * Read a set of input files to determine the bounding box
00271      * of the combined data.
00272      */
00273     static void find_gnuplot_agregate_range( int count,
00274                                              const char* basename,
00275                                              Vector3D& min,
00276                                              Vector3D& max,
00277                                              MsqError& err )
00278     {
00279         // read all input files to determine extents
00280         min = Vector3D( HUGE_VAL, HUGE_VAL, HUGE_VAL );
00281         max = Vector3D( -HUGE_VAL, -HUGE_VAL, -HUGE_VAL );
00282         for( int i = 0; i <= count; ++i )
00283         {
00284             stringstream s;
00285             s << basename << '.' << i << ".gpt";
00286             ifstream infile( s.str().c_str() );
00287             if( !infile )
00288             {
00289                 MSQ_SETERR( err )( s.str(), MsqError::FILE_ACCESS );
00290                 return;
00291             }
00292             double c[3];
00293             while( infile >> c[0] >> c[1] >> c[2] )
00294             {
00295                 for( int j = 0; j < 3; ++j )
00296                 {
00297                     if( c[j] < min[j] ) min[j] = c[j];
00298                     if( c[j] > max[j] ) max[j] = c[j];
00299                 }
00300             }
00301         }
00302     }
00303 
00304     /** Write a GNU Plot script to produce an animation from a
00305      *  sequence of data files
00306      */
00307     void write_gnuplot_animator( int count, const char* basename, MsqError& err )
00308     {
00309         if( count <= 0 ) return;
00310 
00311         const int DELAY  = 10;
00312         const int WIDTH  = 640;
00313         const int HEIGHT = 480;
00314 
00315         // read all input files to determine extents
00316         Vector3D min, max;
00317         find_gnuplot_agregate_range( count, basename, min, max, err );MSQ_ERRRTN( err );
00318 
00319         // chose coordinate plane to plot in
00320         Vector3D range = max - min;
00321         int haxis = 0, vaxis = 1;
00322         if( range[0] < range[1] && range[1] < range[2] )
00323         {
00324             haxis = 1;
00325             vaxis = 2;
00326         }
00327         else if( range[1] < range[2] )
00328         {
00329             vaxis = 2;
00330         }
00331 
00332         // open output file
00333         string base( basename );
00334         ofstream file( ( string( basename ) + ".gnuplot" ).c_str() );
00335         if( !file )
00336         {
00337             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00338             return;
00339         }
00340 
00341         // write header
00342         file << "#!gnuplot" << endl;
00343         file << "#" << endl;
00344         file << "# Mesquite Animation of " << basename << ".0 to " << basename << '.' << count << endl;
00345         file << "#" << endl;
00346 
00347         // write plot settings
00348         file << "set term gif animate transparent opt delay " << DELAY << " size " << WIDTH << "," << HEIGHT << endl;
00349         file << "set xrange [" << min[haxis] - 0.05 * range[haxis] << ":" << max[haxis] + 0.05 * range[haxis] << "]"
00350              << endl;
00351         file << "set yrange [" << min[vaxis] - 0.05 * range[vaxis] << ":" << max[vaxis] + 0.05 * range[vaxis] << "]"
00352              << endl;
00353         file << "set title '" << basename << "'" << endl;
00354         file << "set output '" << basename << ".gif'" << endl;
00355 
00356         // plot data
00357         for( int i = 0; i <= count; ++i )
00358             file << "plot '" << basename << '.' << i << ".gpt'"
00359                  << " using " << haxis + 1 << ":" << vaxis + 1 << " w l" << endl;
00360     }
00361 
00362     static unsigned red( int i, int c )
00363     {
00364         if( i * 4 <= c )
00365             return 255;
00366         else if( i * 4 >= 3 * c )
00367             return 0;
00368         else
00369             return 384 - i * 511 / c;
00370     }
00371 
00372     static unsigned green( int i, int c )
00373     {
00374         if( i * 4 < c )
00375             return i * 510 / c;
00376         else if( i * 4 > 3 * c )
00377             return 1023 - i * 1023 / c;
00378         else
00379             return 255;
00380     }
00381 
00382     static unsigned blue( int i, int c )
00383     {
00384         if( i * 4 <= c )
00385             return 0;
00386         else if( i * 4 >= 3 * c )
00387             return 255;
00388         else
00389             return i * 511 / c - 127;
00390     }
00391 
00392     /** Write a GNU Plot script to produce a single plot from a
00393      *  sequence of data files
00394      */
00395     void write_gnuplot_overlay( int count, const char* basename, MsqError& err )
00396     {
00397         if( count <= 0 ) return;
00398 
00399         const int WIDTH  = 640;
00400         const int HEIGHT = 480;
00401 
00402         // read all input files to determine extents
00403         Vector3D min, max;
00404         find_gnuplot_agregate_range( count, basename, min, max, err );MSQ_ERRRTN( err );
00405 
00406         // chose coordinate plane to plot in
00407         Vector3D range = max - min;
00408         int haxis = 0, vaxis = 1;
00409         if( range[0] < range[1] && range[1] < range[2] )
00410         {
00411             haxis = 1;
00412             vaxis = 2;
00413         }
00414         else if( range[1] < range[2] )
00415         {
00416             vaxis = 2;
00417         }
00418 
00419         // open output file
00420         string base( basename );
00421         FILE* file = fopen( ( string( basename ) + ".gnuplot" ).c_str(), "w" );
00422         if( !file )
00423         {
00424             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00425             return;
00426         }
00427 
00428         // write header
00429         fprintf( file, "#!gnuplot\n" );
00430         fprintf( file, "#\n" );
00431         fprintf( file, "# Mesquite Overlay of %s.0 to %s.%d\n", basename, basename, count );
00432         fprintf( file, "#\n" );
00433 
00434         // write plot settings
00435         fprintf( file, "set term gif size %d,%d\n", WIDTH, HEIGHT );
00436         fprintf( file, "set xrange [%f:%f]\n", min[haxis] - 0.05 * range[haxis], max[haxis] + 0.05 * range[haxis] );
00437         fprintf( file, "set yrange [%f:%f]\n", min[vaxis] - 0.05 * range[vaxis], max[vaxis] + 0.05 * range[vaxis] );
00438         fprintf( file, "set title '%s'\n", basename );
00439         fprintf( file, "set output '%s.gif'\n", basename );
00440 
00441         // plot data
00442         fprintf( file, "plot '%s.0.gpt' using %d:%d w l lc rgb '#%02x%02x%02x' title 't0'", basename, haxis + 1,
00443                  vaxis + 1, red( 0, count ), green( 0, count ), blue( 0, count ) );
00444         for( int i = 1; i <= count; ++i )
00445         {
00446             fprintf( file, ", \\\n     '%s.%d.gpt' using %d:%d w l lc rgb '#%02x%02x%02x' title 't%d'", basename, i,
00447                      haxis + 1, vaxis + 1, red( i, count ), green( i, count ), blue( i, count ), i );
00448         }
00449         fprintf( file, "\n" );
00450 
00451         fclose( file );
00452     }
00453 
00454     void write_stl( Mesh* mesh, const char* filename, MsqError& err )
00455     {
00456         // Open the file
00457         ofstream file( filename );
00458         if( !file )
00459         {
00460             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00461             return;
00462         }
00463 
00464         // Write header
00465         char header[70];
00466         sprintf( header, "Mesquite%d", rand() );
00467         file << "solid " << header << endl;
00468 
00469         MsqVertex coords[3];
00470         std::vector< Mesh::VertexHandle > verts( 3 );
00471         std::vector< size_t > offsets( 2 );
00472 
00473         // Iterate over all elements
00474         size_t count = 0;
00475         std::vector< Mesh::ElementHandle > elems;
00476         std::vector< Mesh::ElementHandle >::iterator iter;
00477         mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
00478         for( iter = elems.begin(); iter != elems.end(); ++iter )
00479         {
00480             // Skip non-triangles
00481             Mesh::ElementHandle elem = *iter;
00482             EntityTopology type;
00483             mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
00484             if( type != TRIANGLE ) continue;
00485             ++count;
00486 
00487             // Get vertex coordinates
00488             mesh->elements_get_attached_vertices( &elem, 1, verts, offsets, err );MSQ_ERRRTN( err );
00489             mesh->vertices_get_coordinates( arrptr( verts ), coords, 3, err );MSQ_ERRRTN( err );
00490 
00491             // Get triagnle normal
00492             Vector3D n = ( coords[0] - coords[1] ) * ( coords[0] - coords[2] );
00493             n.normalize();
00494 
00495             // Write triangle
00496             file << "facet normal " << n.x() << " " << n.y() << " " << n.z() << endl;
00497             file << "outer loop" << endl;
00498             for( unsigned i = 0; i < 3; ++i )
00499                 file << "vertex " << coords[i].x() << " " << coords[i].y() << " " << coords[i].z() << endl;
00500             file << "endloop" << endl;
00501             file << "endfacet" << endl;
00502         }
00503 
00504         file << "endsolid " << header << endl;
00505 
00506         file.close();
00507         if( count == 0 )
00508         {
00509             std::remove( filename );
00510             MSQ_SETERR( err )( "Mesh contains no triangles", MsqError::INVALID_STATE );
00511         }
00512     }
00513 
00514     Projection::Projection( PlanarDomain* domain )
00515     {
00516         Vector3D normal;
00517         domain->vertex_normal_at( 0, normal );
00518         init( normal );
00519     }
00520 
00521     Projection::Projection( const Vector3D& n )
00522     {
00523         init( n );
00524     }
00525 
00526     Projection::Projection( const Vector3D& n, const Vector3D& up )
00527     {
00528         init( n, up );
00529     }
00530 
00531     Projection::Projection( Axis h, Axis v )
00532     {
00533         Vector3D horiz( 0, 0, 0 ), vert( 0, 0, 0 );
00534         horiz[h] = 1.0;
00535         vert[v]  = 1.0;
00536         init( horiz * vert, vert );
00537     }
00538 
00539     void Projection::init( const Vector3D& n )
00540     {
00541         // Choose an "up" direction
00542 
00543         Axis max = X;
00544         for( Axis i = Y; i <= Z; i = (Axis)( i + 1 ) )
00545             if( fabs( n[i] ) > fabs( n[max] ) ) max = i;
00546 
00547         Axis up;
00548         if( max == Y )
00549             up = Z;
00550         else
00551             up = Y;
00552 
00553         // Initialize rotation matrix
00554 
00555         Vector3D up_vect( 0, 0, 0 );
00556         up_vect[up] = 1.0;
00557         init( n, up_vect );
00558     }
00559 
00560     void Projection::init( const Vector3D& n1, const Vector3D& up1 )
00561     {
00562         MsqError err;
00563         const Vector3D n = n1 / n1.length();
00564         const Vector3D u = up1 / up1.length();
00565 
00566         // Rotate for projection
00567         const Vector3D z( 0., 0., 1. );
00568         Vector3D r    = n * z;
00569         double angle  = r.interior_angle( n, z, err );
00570         Matrix3D rot1 = rotation( r, angle );
00571 
00572         // In-plane rotation for up vector
00573         Vector3D pu = u - n * ( n % u );
00574         Vector3D y( 0., 1., 0. );
00575         angle         = z.interior_angle( pu, y, err );
00576         Matrix3D rot2 = rotation( z, angle );
00577 
00578         this->myTransform = rot1 * rot2;
00579     }
00580 
00581     Matrix3D Projection::rotation( const Vector3D& axis, double angle )
00582     {
00583         const double c = cos( angle );
00584         const double s = sin( angle );
00585         const double x = axis.x();
00586         const double y = axis.y();
00587         const double z = axis.z();
00588 
00589         const double xform[9] = { c + x * x * ( 1 - c ),      -z * s + x * y * ( 1 - c ), y * s + x * z * ( 1 - c ),
00590                                   z * s + x * y * ( 1 - c ),  c + y * y * ( 1 - c ),      -x * s + y * z * ( 1 - c ),
00591                                   -y * s + x * z * ( 1 - c ), x * s + y * z * ( 1 - c ),  c + z * z * ( 1 - c ) };
00592         return Matrix3D( xform );
00593     }
00594 
00595     void Projection::project( const Vector3D& p, float& h, float& v )
00596     {
00597         Vector3D pr = myTransform * p;
00598         h           = (float)pr.x();
00599         v           = (float)pr.y();
00600     }
00601 
00602     Transform2D::Transform2D( PatchData* pd, Projection& projection, unsigned width, unsigned height, bool flip )
00603         : myProj( projection ), vertSign( flip ? -1 : 1 )
00604     {
00605         // Get the bounding box of the projected points
00606         float w_max, w_min, h_max, h_min;
00607         w_max = h_max = -std::numeric_limits< float >::max();
00608         w_min = h_min = std::numeric_limits< float >::max();
00609         MsqError err;
00610         const MsqVertex* verts = pd->get_vertex_array( err );
00611         const size_t num_vert  = pd->num_nodes();
00612         for( unsigned i = 0; i < num_vert; ++i )
00613         {
00614             float w, h;
00615             myProj.project( verts[i], w, h );
00616             if( w > w_max ) w_max = w;
00617             if( w < w_min ) w_min = w;
00618             if( h > h_max ) h_max = h;
00619             if( h < h_min ) h_min = h;
00620         }
00621 
00622         // Determine the scale factor
00623         const float w_scale = (float)width / ( w_max - w_min );
00624         const float h_scale = (float)height / ( h_max - h_min );
00625         myScale             = w_scale > h_scale ? h_scale : w_scale;
00626 
00627         // Determine offset
00628         horizOffset = -(int)( myScale * w_min );
00629         vertOffset  = -(int)( myScale * ( flip ? -h_max : h_min ) );
00630 
00631         // Determine bounding box
00632         horizMax = (int)( w_max * myScale ) + horizOffset;
00633         vertMax  = (int)( ( flip ? -h_min : h_max ) * myScale ) + vertOffset;
00634     }
00635 
00636     Transform2D::Transform2D( const Vector3D* verts,
00637                               size_t num_vert,
00638                               Projection& projection,
00639                               unsigned width,
00640                               unsigned height )
00641         : myProj( projection ), vertSign( 1 )
00642     {
00643         // Get the bounding box of the projected points
00644         float w_max, w_min, h_max, h_min;
00645         w_max = h_max = -std::numeric_limits< float >::max();
00646         w_min = h_min = std::numeric_limits< float >::max();
00647         for( unsigned i = 0; i < num_vert; ++i )
00648         {
00649             float w, h;
00650             myProj.project( verts[i], w, h );
00651             if( w > w_max ) w_max = w;
00652             if( w < w_min ) w_min = w;
00653             if( h > h_max ) h_max = h;
00654             if( h < h_min ) h_min = h;
00655         }
00656 
00657         // Determine the scale factor
00658         const float w_scale = (float)width / ( w_max - w_min );
00659         const float h_scale = (float)height / ( h_max - h_min );
00660         myScale             = w_scale > h_scale ? h_scale : w_scale;
00661 
00662         // Determine offset
00663         horizOffset = -(int)( myScale * w_min );
00664         vertOffset  = -(int)( myScale * h_min );
00665 
00666         // Determine bounding box
00667         horizMax = (int)( w_max * myScale ) + horizOffset;
00668         vertMax  = (int)( h_max * myScale ) + vertOffset;
00669     }
00670 
00671     void Transform2D::transform( const Vector3D& coords, int& horizontal, int& vertical ) const
00672     {
00673         float horiz, vert;
00674         myProj.project( coords, horiz, vert );
00675         horizontal = (int)( myScale * horiz ) + horizOffset;
00676         vertical   = vertSign * (int)( myScale * vert ) + vertOffset;
00677     }
00678 
00679     /** Write quadratic edge shape in PostScript format.
00680      *
00681      * Given the three points composing a quadratic mesh edge,
00682      * write the cubic Bezier curve of the same shape in
00683      * PostScript format.  The formulas for P1 and P2
00684      * at the start of this function will result in the cubic
00685      * terms of the Bezier curve dropping out, leaving the
00686      * quadratic curve matching the edge shape function as
00687      * described in Section 3.6 of Hughes.  (If you're attempting
00688      * to verify this, don't forget to adjust for the different
00689      * parameter ranges: \f$ \xi = 2 t - 1 \f$).
00690      */
00691     static void write_eps_quadratic_edge( ostream& s, Transform2D& xform, Vector3D start, Vector3D mid, Vector3D end )
00692     {
00693         Vector3D P1 = 1. / 3 * ( 4 * mid - end );
00694         Vector3D P2 = 1. / 3 * ( 4 * mid - start );
00695 
00696         int x, y;
00697         xform.transform( start, x, y );
00698         s << x << ' ' << y << " moveto" << endl;
00699         xform.transform( P1, x, y );
00700         s << x << ' ' << y << ' ';
00701         xform.transform( P2, x, y );
00702         s << x << ' ' << y << ' ';
00703         xform.transform( end, x, y );
00704         s << x << ' ' << y << " curveto" << endl;
00705     }
00706 
00707     void write_eps( Mesh* mesh, const char* filename, Projection proj, MsqError& err, int width, int height )
00708     {
00709         // Get a global patch
00710         PatchData pd;
00711         pd.set_mesh( mesh );
00712         pd.fill_global_patch( err );MSQ_ERRRTN( err );
00713 
00714         Transform2D transf( &pd, proj, width, height, false );
00715 
00716         // Open the file
00717         ofstream s( filename );
00718         if( !s )
00719         {
00720             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00721             return;
00722         }
00723 
00724         // Write header
00725         s << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
00726         s << "%%Creator: Mesquite" << endl;
00727         s << "%%Title: Mesquite " << endl;
00728         s << "%%DocumentData: Clean7Bit" << endl;
00729         s << "%%Origin: 0 0" << endl;
00730         s << "%%BoundingBox: 0 0 " << transf.max_horizontal() << ' ' << transf.max_vertical() << endl;
00731         s << "%%Pages: 1" << endl;
00732 
00733         s << "%%BeginProlog" << endl;
00734         s << "save" << endl;
00735         s << "countdictstack" << endl;
00736         s << "mark" << endl;
00737         s << "newpath" << endl;
00738         s << "/showpage {} def" << endl;
00739         s << "/setpagedevice {pop} def" << endl;
00740         s << "%%EndProlog" << endl;
00741 
00742         s << "%%Page: 1 1" << endl;
00743         s << "1 setlinewidth" << endl;
00744         s << "0.0 setgray" << endl;
00745 
00746         // Write mesh edges
00747         EdgeIterator iter( &pd, err );MSQ_ERRRTN( err );
00748         while( !iter.is_at_end() )
00749         {
00750             int s_w, s_h, e_w, e_h;
00751             transf.transform( iter.start(), s_w, s_h );
00752             transf.transform( iter.end(), e_w, e_h );
00753 
00754             s << "newpath" << endl;
00755             s << s_w << ' ' << s_h << " moveto" << endl;
00756 
00757             if( !iter.mid() )
00758             {
00759                 s << e_w << ' ' << e_h << " lineto" << endl;
00760             }
00761             else
00762             {
00763                 write_eps_quadratic_edge( s, transf, iter.start(), *iter.mid(), iter.end() );
00764                 // draw rings at mid-edge node location
00765                 // transf.transform( *(iter.mid()), w1, h1 );
00766                 // s << w1+2 << ' ' << h1 <<  " moveto"            << endl;
00767                 // s << w1 << ' ' << h1 <<  " 2 0 360 arc"         << endl;
00768             }
00769             s << "stroke" << endl;
00770 
00771             iter.step( err );MSQ_ERRRTN( err );
00772         }
00773 
00774         // Write footer
00775         s << "%%Trailer" << endl;
00776         s << "cleartomark" << endl;
00777         s << "countdictstack" << endl;
00778         s << "exch sub { end } repeat" << endl;
00779         s << "restore" << endl;
00780         s << "%%EOF" << endl;
00781     }
00782 
00783     /** Quadratic triangle shape function for use in write_eps_triangle */
00784     static double tN0( double r, double s )
00785     {
00786         double t = 1 - r - s;
00787         return t * ( 2 * t - 1 );
00788     }
00789     static double tN1( double r, double )
00790     {
00791         return r * ( 2 * r - 1 );
00792     }
00793     static double tN2( double, double s )
00794     {
00795         return s * ( 2 * s - 1 );
00796     }
00797     static double tN3( double r, double s )
00798     {
00799         double t = 1 - r - s;
00800         return 4 * r * t;
00801     }
00802     static double tN4( double r, double s )
00803     {
00804         return 4 * r * s;
00805     }
00806     static double tN5( double r, double s )
00807     {
00808         double t = 1 - r - s;
00809         return 4 * s * t;
00810     }
00811     /** Quadratic triangle shape function for use in write_eps_triangle */
00812     static Vector3D quad_tri_pt( double r, double s, const Vector3D* coords )
00813     {
00814         Vector3D result = tN0( r, s ) * coords[0];
00815         result += tN1( r, s ) * coords[1];
00816         result += tN2( r, s ) * coords[2];
00817         result += tN3( r, s ) * coords[3];
00818         result += tN4( r, s ) * coords[4];
00819         result += tN5( r, s ) * coords[5];
00820         return result;
00821     }
00822 
00823     void write_eps_triangle( Mesh* mesh,
00824                              Mesh::ElementHandle elem,
00825                              const char* filename,
00826                              bool draw_iso_lines,
00827                              bool draw_nodes,
00828                              MsqError& err,
00829                              int width,
00830                              int height )
00831     {
00832         // Get triangle vertices
00833         MsqVertex coords[6];
00834         EntityTopology type;
00835         mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
00836         if( type != TRIANGLE )
00837         {
00838             MSQ_SETERR( err )( "Invalid element type", MsqError::UNSUPPORTED_ELEMENT );
00839             return;
00840         }
00841         std::vector< Mesh::VertexHandle > verts;
00842         std::vector< size_t > junk;
00843         mesh->elements_get_attached_vertices( &elem, 1, verts, junk, err );MSQ_ERRRTN( err );
00844         if( verts.size() != 3 && verts.size() != 6 )
00845         {
00846             MSQ_SETERR( err )( "Invalid element type", MsqError::UNSUPPORTED_ELEMENT );
00847             return;
00848         }
00849         mesh->vertices_get_coordinates( arrptr( verts ), coords, verts.size(), err );MSQ_ERRRTN( err );
00850 
00851         Vector3D coords2[6];
00852         std::copy( coords, coords + verts.size(), coords2 );
00853 
00854         std::vector< bool > fixed( verts.size(), false );
00855         if( draw_nodes )
00856         {
00857             mesh->vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err );MSQ_ERRRTN( err );
00858         }
00859         write_eps_triangle( coords2, verts.size(), filename, draw_iso_lines, draw_nodes, err, fixed, width, height );
00860     }
00861 
00862     void write_eps_triangle( const Vector3D* coords,
00863                              size_t num_vtx,
00864                              const char* filename,
00865                              bool draw_iso_lines,
00866                              bool draw_nodes,
00867                              MsqError& err,
00868                              const std::vector< bool >& fixed,
00869                              int width,
00870                              int height )
00871     {
00872         const int PT_RAD        = 3;           // radius of circles for drawing nodes, in points
00873         const int PAD           = PT_RAD + 2;  // margin in points
00874         const double EDGE_GRAY  = 0.0;         // color for triangle edges, 0.0 => black
00875         const double ISO_GRAY   = 0.7;         // color for parameter iso-lines
00876         const double NODE_GRAY  = 0.0;         // color for node circle
00877         const double FIXED_GRAY = 1.0;         // color to fill fixed nodes with, 1.0 => white
00878         const double FREE_GRAY  = 0.0;         // color to fill free nodes with
00879 
00880         Projection proj( X, Y );
00881         Transform2D transf( coords, num_vtx, proj, width, height );
00882 
00883         // Open the file
00884         ofstream str( filename );
00885         if( !str )
00886         {
00887             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00888             return;
00889         }
00890 
00891         // Write header
00892         str << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
00893         str << "%%Creator: Mesquite" << endl;
00894         str << "%%Title: Mesquite " << endl;
00895         str << "%%DocumentData: Clean7Bit" << endl;
00896         str << "%%Origin: 0 0" << endl;
00897         str << "%%BoundingBox: " << -PAD << ' ' << -PAD << ' ' << transf.max_horizontal() + PAD << ' '
00898             << transf.max_vertical() + PAD << endl;
00899         str << "%%Pages: 1" << endl;
00900 
00901         str << "%%BeginProlog" << endl;
00902         str << "save" << endl;
00903         str << "countdictstack" << endl;
00904         str << "mark" << endl;
00905         str << "newpath" << endl;
00906         str << "/showpage {} def" << endl;
00907         str << "/setpagedevice {pop} def" << endl;
00908         str << "%%EndProlog" << endl;
00909 
00910         str << "%%Page: 1 1" << endl;
00911         str << "1 setlinewidth" << endl;
00912         str << EDGE_GRAY << " setgray" << endl;
00913 
00914         const double h = 0.5, t = 1.0 / 3.0, w = 2.0 / 3.0, s = 1. / 6, f = 5. / 6;
00915         const int NUM_ISO                      = 15;
00916         const double iso_params[NUM_ISO][2][2] = {
00917             { { h, 0 }, { h, h } },  // r = 1/2
00918             { { t, 0 }, { t, w } },  // r = 1/3
00919             { { w, 0 }, { w, t } },  // r = 2/3
00920             { { s, 0 }, { s, f } },  // r = 1/6
00921             { { f, 0 }, { f, s } },  // r = 5/6
00922             { { 0, h }, { h, h } },  // s = 1/2
00923             { { 0, t }, { w, t } },  // s = 1/3
00924             { { 0, w }, { t, w } },  // s = 2/3
00925             { { 0, s }, { f, s } },  // s = 1/6
00926             { { 0, f }, { s, f } },  // s = 5/6
00927             { { 0, h }, { h, 0 } },  // t = 1 - r - s = 1/2
00928             { { 0, w }, { w, 0 } },  // t = 1 - r - s = 1/3
00929             { { 0, t }, { t, 0 } },  // t = 1 - r - s = 2/3
00930             { { 0, f }, { f, 0 } },  // t = 1 - r - s = 1/6
00931             { { 0, s }, { s, 0 } }   // t = 1 - r - s = 5/6
00932         };
00933 
00934         if( num_vtx == 3 )
00935         {
00936             int x[3], y[3];
00937             for( size_t i = 0; i < 3; ++i )
00938                 transf.transform( coords[i], x[i], y[i] );
00939 
00940             str << "newpath" << endl;
00941             str << x[0] << ' ' << y[0] << " moveto" << endl;
00942             str << x[1] << ' ' << y[1] << " lineto" << endl;
00943             str << x[2] << ' ' << y[2] << " lineto" << endl;
00944             str << x[0] << ' ' << y[0] << " lineto" << endl;
00945             str << "stroke" << endl;
00946 
00947             if( draw_iso_lines )
00948             {
00949                 str << ISO_GRAY << " setgray" << endl;
00950                 str << "newpath" << endl;
00951                 for( int i = 0; i < NUM_ISO; ++i )
00952                 {
00953                     double R[2]   = { iso_params[i][0][0], iso_params[i][1][0] };
00954                     double S[2]   = { iso_params[i][0][1], iso_params[i][1][1] };
00955                     double T[2]   = { 1 - R[0] - S[0], 1 - R[1] - S[1] };
00956                     Vector3D p[2] = { T[0] * coords[0] + R[0] * coords[1] + S[0] * coords[2],
00957                                       T[1] * coords[0] + R[1] * coords[1] + S[1] * coords[2] };
00958                     transf.transform( p[0], x[0], y[0] );
00959                     transf.transform( p[1], x[1], y[1] );
00960                     str << x[0] << ' ' << y[0] << " moveto" << endl;
00961                     str << x[1] << ' ' << y[1] << " lineto" << endl;
00962                 }
00963                 str << "    stroke" << endl;
00964             }
00965         }
00966         else if( num_vtx == 6 )
00967         {
00968             str << "newpath" << endl;
00969             write_eps_quadratic_edge( str, transf, coords[0], coords[3], coords[1] );
00970             write_eps_quadratic_edge( str, transf, coords[1], coords[4], coords[2] );
00971             write_eps_quadratic_edge( str, transf, coords[2], coords[5], coords[0] );
00972             str << "stroke" << endl;
00973 
00974             if( draw_iso_lines )
00975             {
00976                 str << ISO_GRAY << " setgray" << endl;
00977                 str << "newpath" << endl;
00978                 for( int i = 0; i < NUM_ISO; ++i )
00979                 {
00980                     double R[3]   = { iso_params[i][0][0], 0, iso_params[i][1][0] };
00981                     double S[3]   = { iso_params[i][0][1], 0, iso_params[i][1][1] };
00982                     R[1]          = 0.5 * ( R[0] + R[2] );
00983                     S[1]          = 0.5 * ( S[0] + S[2] );
00984                     Vector3D p[3] = { quad_tri_pt( R[0], S[0], coords ), quad_tri_pt( R[1], S[1], coords ),
00985                                       quad_tri_pt( R[2], S[2], coords ) };
00986                     write_eps_quadratic_edge( str, transf, p[0], p[1], p[2] );
00987                 }
00988                 str << "    stroke" << endl;
00989             }
00990         }
00991 
00992         if( draw_nodes )
00993         {
00994             for( size_t i = 0; i < num_vtx; ++i )
00995             {
00996                 int iw, ih;
00997                 // fill interior with either white or black depending
00998                 // on whether or not the vertex is fixed.
00999                 if( fixed[i] )
01000                     str << FIXED_GRAY << " setgray" << endl;
01001                 else
01002                     str << FREE_GRAY << " setgray" << endl;
01003                 transf.transform( coords[i], iw, ih );
01004                 str << w + PT_RAD << ' ' << ih << " moveto" << endl;
01005                 str << iw << ' ' << ih << ' ' << PT_RAD << " 0 360 arc" << endl;
01006                 str << "closepath fill" << endl;
01007                 str << NODE_GRAY << " setgray" << endl;
01008                 str << "newpath" << endl;
01009                 str << iw << ' ' << ih << ' ' << PT_RAD << " 0 360 arc" << endl;
01010                 str << "stroke" << endl;
01011             }
01012         }
01013 
01014         // Write footer
01015         str << "%%Trailer" << endl;
01016         str << "cleartomark" << endl;
01017         str << "countdictstack" << endl;
01018         str << "exch sub { end } repeat" << endl;
01019         str << "restore" << endl;
01020         str << "%%EOF" << endl;
01021     }
01022 
01023     void write_svg( Mesh* mesh, const char* filename, Projection proj, MsqError& err )
01024     {
01025         // Get a global patch
01026         PatchData pd;
01027         pd.set_mesh( mesh );
01028         pd.fill_global_patch( err );MSQ_ERRRTN( err );
01029 
01030         Transform2D transf( &pd, proj, 400, 400, true );
01031 
01032         // Open the file
01033         ofstream file( filename );
01034         if( !file )
01035         {
01036             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
01037             return;
01038         }
01039 
01040         // Write header
01041         file << "<?xml version=\"1.0\" standalone=\"no\"?>" << endl;
01042         file << "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" "
01043              << "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << endl;
01044         file << endl;
01045         file << "<svg width=\"" << transf.max_horizontal() << "\" height=\"" << transf.max_vertical()
01046              << "\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" << endl;
01047 
01048         // Write mesh edges
01049         EdgeIterator iter( &pd, err );MSQ_ERRRTN( err );
01050         while( !iter.is_at_end() )
01051         {
01052             int s_w, s_h, e_w, e_h;
01053             transf.transform( iter.start(), s_w, s_h );
01054             transf.transform( iter.end(), e_w, e_h );
01055 
01056             file << "<line "
01057                  << "x1=\"" << s_w << "\" "
01058                  << "y1=\"" << s_h << "\" "
01059                  << "x2=\"" << e_w << "\" "
01060                  << "y2=\"" << e_h << "\" "
01061                  << " style=\"stroke:rgb(99,99,99);stroke-width:2\""
01062                  << "/>" << endl;
01063 
01064             iter.step( err );MSQ_ERRRTN( err );
01065         }
01066 
01067         // Write footer
01068         file << "</svg>" << endl;
01069     }
01070 
01071 }  // namespace MeshWriter
01072 
01073 }  // namespace MBMesquite
01074 
01075 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines