MOAB: Mesh Oriented datABase  (version 5.2.1)
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, const char* basename, Vector3D& min, Vector3D& max,
00274                                              MsqError& err )
00275     {
00276         // read all input files to determine extents
00277         min = Vector3D( HUGE_VAL, HUGE_VAL, HUGE_VAL );
00278         max = Vector3D( -HUGE_VAL, -HUGE_VAL, -HUGE_VAL );
00279         for( int i = 0; i <= count; ++i )
00280         {
00281             stringstream s;
00282             s << basename << '.' << i << ".gpt";
00283             ifstream infile( s.str().c_str() );
00284             if( !infile )
00285             {
00286                 MSQ_SETERR( err )( s.str(), MsqError::FILE_ACCESS );
00287                 return;
00288             }
00289             double c[3];
00290             while( infile >> c[0] >> c[1] >> c[2] )
00291             {
00292                 for( int j = 0; j < 3; ++j )
00293                 {
00294                     if( c[j] < min[j] ) min[j] = c[j];
00295                     if( c[j] > max[j] ) max[j] = c[j];
00296                 }
00297             }
00298         }
00299     }
00300 
00301     /** Write a GNU Plot script to produce an animation from a
00302      *  sequence of data files
00303      */
00304     void write_gnuplot_animator( int count, const char* basename, MsqError& err )
00305     {
00306         if( count <= 0 ) return;
00307 
00308         const int DELAY  = 10;
00309         const int WIDTH  = 640;
00310         const int HEIGHT = 480;
00311 
00312         // read all input files to determine extents
00313         Vector3D min, max;
00314         find_gnuplot_agregate_range( count, basename, min, max, err );MSQ_ERRRTN( err );
00315 
00316         // chose coordinate plane to plot in
00317         Vector3D range = max - min;
00318         int haxis = 0, vaxis = 1;
00319         if( range[0] < range[1] && range[1] < range[2] )
00320         {
00321             haxis = 1;
00322             vaxis = 2;
00323         }
00324         else if( range[1] < range[2] )
00325         {
00326             vaxis = 2;
00327         }
00328 
00329         // open output file
00330         string base( basename );
00331         ofstream file( ( string( basename ) + ".gnuplot" ).c_str() );
00332         if( !file )
00333         {
00334             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00335             return;
00336         }
00337 
00338         // write header
00339         file << "#!gnuplot" << endl;
00340         file << "#" << endl;
00341         file << "# Mesquite Animation of " << basename << ".0 to " << basename << '.' << count << endl;
00342         file << "#" << endl;
00343 
00344         // write plot settings
00345         file << "set term gif animate transparent opt delay " << DELAY << " size " << WIDTH << "," << HEIGHT << endl;
00346         file << "set xrange [" << min[haxis] - 0.05 * range[haxis] << ":" << max[haxis] + 0.05 * range[haxis] << "]"
00347              << endl;
00348         file << "set yrange [" << min[vaxis] - 0.05 * range[vaxis] << ":" << max[vaxis] + 0.05 * range[vaxis] << "]"
00349              << endl;
00350         file << "set title '" << basename << "'" << endl;
00351         file << "set output '" << basename << ".gif'" << endl;
00352 
00353         // plot data
00354         for( int i = 0; i <= count; ++i )
00355             file << "plot '" << basename << '.' << i << ".gpt'"
00356                  << " using " << haxis + 1 << ":" << vaxis + 1 << " w l" << endl;
00357     }
00358 
00359     static unsigned red( int i, int c )
00360     {
00361         if( i * 4 <= c )
00362             return 255;
00363         else if( i * 4 >= 3 * c )
00364             return 0;
00365         else
00366             return 384 - i * 511 / c;
00367     }
00368 
00369     static unsigned green( int i, int c )
00370     {
00371         if( i * 4 < c )
00372             return i * 510 / c;
00373         else if( i * 4 > 3 * c )
00374             return 1023 - i * 1023 / c;
00375         else
00376             return 255;
00377     }
00378 
00379     static unsigned blue( int i, int c )
00380     {
00381         if( i * 4 <= c )
00382             return 0;
00383         else if( i * 4 >= 3 * c )
00384             return 255;
00385         else
00386             return i * 511 / c - 127;
00387     }
00388 
00389     /** Write a GNU Plot script to produce a single plot from a
00390      *  sequence of data files
00391      */
00392     void write_gnuplot_overlay( int count, const char* basename, MsqError& err )
00393     {
00394         if( count <= 0 ) return;
00395 
00396         const int WIDTH  = 640;
00397         const int HEIGHT = 480;
00398 
00399         // read all input files to determine extents
00400         Vector3D min, max;
00401         find_gnuplot_agregate_range( count, basename, min, max, err );MSQ_ERRRTN( err );
00402 
00403         // chose coordinate plane to plot in
00404         Vector3D range = max - min;
00405         int haxis = 0, vaxis = 1;
00406         if( range[0] < range[1] && range[1] < range[2] )
00407         {
00408             haxis = 1;
00409             vaxis = 2;
00410         }
00411         else if( range[1] < range[2] )
00412         {
00413             vaxis = 2;
00414         }
00415 
00416         // open output file
00417         string base( basename );
00418         FILE* file = fopen( ( string( basename ) + ".gnuplot" ).c_str(), "w" );
00419         if( !file )
00420         {
00421             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00422             return;
00423         }
00424 
00425         // write header
00426         fprintf( file, "#!gnuplot\n" );
00427         fprintf( file, "#\n" );
00428         fprintf( file, "# Mesquite Overlay of %s.0 to %s.%d\n", basename, basename, count );
00429         fprintf( file, "#\n" );
00430 
00431         // write plot settings
00432         fprintf( file, "set term gif size %d,%d\n", WIDTH, HEIGHT );
00433         fprintf( file, "set xrange [%f:%f]\n", min[haxis] - 0.05 * range[haxis], max[haxis] + 0.05 * range[haxis] );
00434         fprintf( file, "set yrange [%f:%f]\n", min[vaxis] - 0.05 * range[vaxis], max[vaxis] + 0.05 * range[vaxis] );
00435         fprintf( file, "set title '%s'\n", basename );
00436         fprintf( file, "set output '%s.gif'\n", basename );
00437 
00438         // plot data
00439         fprintf( file, "plot '%s.0.gpt' using %d:%d w l lc rgb '#%02x%02x%02x' title 't0'", basename, haxis + 1,
00440                  vaxis + 1, red( 0, count ), green( 0, count ), blue( 0, count ) );
00441         for( int i = 1; i <= count; ++i )
00442         {
00443             fprintf( file, ", \\\n     '%s.%d.gpt' using %d:%d w l lc rgb '#%02x%02x%02x' title 't%d'", basename, i,
00444                      haxis + 1, vaxis + 1, red( i, count ), green( i, count ), blue( i, count ), i );
00445         }
00446         fprintf( file, "\n" );
00447 
00448         fclose( file );
00449     }
00450 
00451     void write_stl( Mesh* mesh, const char* filename, MsqError& err )
00452     {
00453         // Open the file
00454         ofstream file( filename );
00455         if( !file )
00456         {
00457             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00458             return;
00459         }
00460 
00461         // Write header
00462         char header[70];
00463         sprintf( header, "Mesquite%d", rand() );
00464         file << "solid " << header << endl;
00465 
00466         MsqVertex coords[3];
00467         std::vector< Mesh::VertexHandle > verts( 3 );
00468         std::vector< size_t > offsets( 2 );
00469 
00470         // Iterate over all elements
00471         size_t count = 0;
00472         std::vector< Mesh::ElementHandle > elems;
00473         std::vector< Mesh::ElementHandle >::iterator iter;
00474         mesh->get_all_elements( elems, err );MSQ_ERRRTN( err );
00475         for( iter = elems.begin(); iter != elems.end(); ++iter )
00476         {
00477             // Skip non-triangles
00478             Mesh::ElementHandle elem = *iter;
00479             EntityTopology type;
00480             mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
00481             if( type != TRIANGLE ) continue;
00482             ++count;
00483 
00484             // Get vertex coordinates
00485             mesh->elements_get_attached_vertices( &elem, 1, verts, offsets, err );MSQ_ERRRTN( err );
00486             mesh->vertices_get_coordinates( arrptr( verts ), coords, 3, err );MSQ_ERRRTN( err );
00487 
00488             // Get triagnle normal
00489             Vector3D n = ( coords[0] - coords[1] ) * ( coords[0] - coords[2] );
00490             n.normalize();
00491 
00492             // Write triangle
00493             file << "facet normal " << n.x() << " " << n.y() << " " << n.z() << endl;
00494             file << "outer loop" << endl;
00495             for( unsigned i = 0; i < 3; ++i )
00496                 file << "vertex " << coords[i].x() << " " << coords[i].y() << " " << coords[i].z() << endl;
00497             file << "endloop" << endl;
00498             file << "endfacet" << endl;
00499         }
00500 
00501         file << "endsolid " << header << endl;
00502 
00503         file.close();
00504         if( count == 0 )
00505         {
00506             std::remove( filename );
00507             MSQ_SETERR( err )( "Mesh contains no triangles", MsqError::INVALID_STATE );
00508         }
00509     }
00510 
00511     Projection::Projection( PlanarDomain* domain )
00512     {
00513         Vector3D normal;
00514         domain->vertex_normal_at( 0, normal );
00515         init( normal );
00516     }
00517 
00518     Projection::Projection( const Vector3D& n )
00519     {
00520         init( n );
00521     }
00522 
00523     Projection::Projection( const Vector3D& n, const Vector3D& up )
00524     {
00525         init( n, up );
00526     }
00527 
00528     Projection::Projection( Axis h, Axis v )
00529     {
00530         Vector3D horiz( 0, 0, 0 ), vert( 0, 0, 0 );
00531         horiz[h] = 1.0;
00532         vert[v]  = 1.0;
00533         init( horiz * vert, vert );
00534     }
00535 
00536     void Projection::init( const Vector3D& n )
00537     {
00538         // Choose an "up" direction
00539 
00540         Axis max = X;
00541         for( Axis i = Y; i <= Z; i = ( Axis )( i + 1 ) )
00542             if( fabs( n[i] ) > fabs( n[max] ) ) max = i;
00543 
00544         Axis up;
00545         if( max == Y )
00546             up = Z;
00547         else
00548             up = Y;
00549 
00550         // Initialize rotation matrix
00551 
00552         Vector3D up_vect( 0, 0, 0 );
00553         up_vect[up] = 1.0;
00554         init( n, up_vect );
00555     }
00556 
00557     void Projection::init( const Vector3D& n1, const Vector3D& up1 )
00558     {
00559         MsqError err;
00560         const Vector3D n = n1 / n1.length();
00561         const Vector3D u = up1 / up1.length();
00562 
00563         // Rotate for projection
00564         const Vector3D z( 0., 0., 1. );
00565         Vector3D r    = n * z;
00566         double angle  = r.interior_angle( n, z, err );
00567         Matrix3D rot1 = rotation( r, angle );
00568 
00569         // In-plane rotation for up vector
00570         Vector3D pu = u - n * ( n % u );
00571         Vector3D y( 0., 1., 0. );
00572         angle         = z.interior_angle( pu, y, err );
00573         Matrix3D rot2 = rotation( z, angle );
00574 
00575         this->myTransform = rot1 * rot2;
00576     }
00577 
00578     Matrix3D Projection::rotation( const Vector3D& axis, double angle )
00579     {
00580         const double c = cos( angle );
00581         const double s = sin( angle );
00582         const double x = axis.x();
00583         const double y = axis.y();
00584         const double z = axis.z();
00585 
00586         const double xform[9] = { c + x * x * ( 1 - c ),      -z * s + x * y * ( 1 - c ), y * s + x * z * ( 1 - c ),
00587                                   z * s + x * y * ( 1 - c ),  c + y * y * ( 1 - c ),      -x * s + y * z * ( 1 - c ),
00588                                   -y * s + x * z * ( 1 - c ), x * s + y * z * ( 1 - c ),  c + z * z * ( 1 - c ) };
00589         return Matrix3D( xform );
00590     }
00591 
00592     void Projection::project( const Vector3D& p, float& h, float& v )
00593     {
00594         Vector3D pr = myTransform * p;
00595         h           = (float)pr.x();
00596         v           = (float)pr.y();
00597     }
00598 
00599     Transform2D::Transform2D( PatchData* pd, Projection& projection, unsigned width, unsigned height, bool flip )
00600         : myProj( projection ), vertSign( flip ? -1 : 1 )
00601     {
00602         // Get the bounding box of the projected points
00603         float w_max, w_min, h_max, h_min;
00604         w_max = h_max = -std::numeric_limits< float >::max();
00605         w_min = h_min = std::numeric_limits< float >::max();
00606         MsqError err;
00607         const MsqVertex* verts = pd->get_vertex_array( err );
00608         const size_t num_vert  = pd->num_nodes();
00609         for( unsigned i = 0; i < num_vert; ++i )
00610         {
00611             float w, h;
00612             myProj.project( verts[i], w, h );
00613             if( w > w_max ) w_max = w;
00614             if( w < w_min ) w_min = w;
00615             if( h > h_max ) h_max = h;
00616             if( h < h_min ) h_min = h;
00617         }
00618 
00619         // Determine the scale factor
00620         const float w_scale = (float)width / ( w_max - w_min );
00621         const float h_scale = (float)height / ( h_max - h_min );
00622         myScale             = w_scale > h_scale ? h_scale : w_scale;
00623 
00624         // Determine offset
00625         horizOffset = -(int)( myScale * w_min );
00626         vertOffset  = -(int)( myScale * ( flip ? -h_max : h_min ) );
00627 
00628         // Determine bounding box
00629         horizMax = (int)( w_max * myScale ) + horizOffset;
00630         vertMax  = (int)( ( flip ? -h_min : h_max ) * myScale ) + vertOffset;
00631     }
00632 
00633     Transform2D::Transform2D( const Vector3D* verts, size_t num_vert, Projection& projection, unsigned width,
00634                               unsigned height )
00635         : myProj( projection ), vertSign( 1 )
00636     {
00637         // Get the bounding box of the projected points
00638         float w_max, w_min, h_max, h_min;
00639         w_max = h_max = -std::numeric_limits< float >::max();
00640         w_min = h_min = std::numeric_limits< float >::max();
00641         for( unsigned i = 0; i < num_vert; ++i )
00642         {
00643             float w, h;
00644             myProj.project( verts[i], w, h );
00645             if( w > w_max ) w_max = w;
00646             if( w < w_min ) w_min = w;
00647             if( h > h_max ) h_max = h;
00648             if( h < h_min ) h_min = h;
00649         }
00650 
00651         // Determine the scale factor
00652         const float w_scale = (float)width / ( w_max - w_min );
00653         const float h_scale = (float)height / ( h_max - h_min );
00654         myScale             = w_scale > h_scale ? h_scale : w_scale;
00655 
00656         // Determine offset
00657         horizOffset = -(int)( myScale * w_min );
00658         vertOffset  = -(int)( myScale * h_min );
00659 
00660         // Determine bounding box
00661         horizMax = (int)( w_max * myScale ) + horizOffset;
00662         vertMax  = (int)( h_max * myScale ) + vertOffset;
00663     }
00664 
00665     void Transform2D::transform( const Vector3D& coords, int& horizontal, int& vertical ) const
00666     {
00667         float horiz, vert;
00668         myProj.project( coords, horiz, vert );
00669         horizontal = (int)( myScale * horiz ) + horizOffset;
00670         vertical   = vertSign * (int)( myScale * vert ) + vertOffset;
00671     }
00672 
00673     /** Write quadratic edge shape in PostScript format.
00674      *
00675      * Given the three points composing a quadratic mesh edge,
00676      * write the cubic Bezier curve of the same shape in
00677      * PostScript format.  The formulas for P1 and P2
00678      * at the start of this function will result in the cubic
00679      * terms of the Bezier curve dropping out, leaving the
00680      * quadratic curve matching the edge shape function as
00681      * described in Section 3.6 of Hughes.  (If you're attempting
00682      * to verify this, don't forget to adjust for the different
00683      * parameter ranges: \f$ \xi = 2 t - 1 \f$).
00684      */
00685     static void write_eps_quadratic_edge( ostream& s, Transform2D& xform, Vector3D start, Vector3D mid, Vector3D end )
00686     {
00687         Vector3D P1 = 1. / 3 * ( 4 * mid - end );
00688         Vector3D P2 = 1. / 3 * ( 4 * mid - start );
00689 
00690         int x, y;
00691         xform.transform( start, x, y );
00692         s << x << ' ' << y << " moveto" << endl;
00693         xform.transform( P1, x, y );
00694         s << x << ' ' << y << ' ';
00695         xform.transform( P2, x, y );
00696         s << x << ' ' << y << ' ';
00697         xform.transform( end, x, y );
00698         s << x << ' ' << y << " curveto" << endl;
00699     }
00700 
00701     void write_eps( Mesh* mesh, const char* filename, Projection proj, MsqError& err, int width, int height )
00702     {
00703         // Get a global patch
00704         PatchData pd;
00705         pd.set_mesh( mesh );
00706         pd.fill_global_patch( err );MSQ_ERRRTN( err );
00707 
00708         Transform2D transf( &pd, proj, width, height, false );
00709 
00710         // Open the file
00711         ofstream s( filename );
00712         if( !s )
00713         {
00714             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00715             return;
00716         }
00717 
00718         // Write header
00719         s << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
00720         s << "%%Creator: Mesquite" << endl;
00721         s << "%%Title: Mesquite " << endl;
00722         s << "%%DocumentData: Clean7Bit" << endl;
00723         s << "%%Origin: 0 0" << endl;
00724         s << "%%BoundingBox: 0 0 " << transf.max_horizontal() << ' ' << transf.max_vertical() << endl;
00725         s << "%%Pages: 1" << endl;
00726 
00727         s << "%%BeginProlog" << endl;
00728         s << "save" << endl;
00729         s << "countdictstack" << endl;
00730         s << "mark" << endl;
00731         s << "newpath" << endl;
00732         s << "/showpage {} def" << endl;
00733         s << "/setpagedevice {pop} def" << endl;
00734         s << "%%EndProlog" << endl;
00735 
00736         s << "%%Page: 1 1" << endl;
00737         s << "1 setlinewidth" << endl;
00738         s << "0.0 setgray" << endl;
00739 
00740         // Write mesh edges
00741         EdgeIterator iter( &pd, err );MSQ_ERRRTN( err );
00742         while( !iter.is_at_end() )
00743         {
00744             int s_w, s_h, e_w, e_h;
00745             transf.transform( iter.start(), s_w, s_h );
00746             transf.transform( iter.end(), e_w, e_h );
00747 
00748             s << "newpath" << endl;
00749             s << s_w << ' ' << s_h << " moveto" << endl;
00750 
00751             if( !iter.mid() ) { s << e_w << ' ' << e_h << " lineto" << endl; }
00752             else
00753             {
00754                 write_eps_quadratic_edge( s, transf, iter.start(), *iter.mid(), iter.end() );
00755                 // draw rings at mid-edge node location
00756                 // transf.transform( *(iter.mid()), w1, h1 );
00757                 // s << w1+2 << ' ' << h1 <<  " moveto"            << endl;
00758                 // s << w1 << ' ' << h1 <<  " 2 0 360 arc"         << endl;
00759             }
00760             s << "stroke" << endl;
00761 
00762             iter.step( err );MSQ_ERRRTN( err );
00763         }
00764 
00765         // Write footer
00766         s << "%%Trailer" << endl;
00767         s << "cleartomark" << endl;
00768         s << "countdictstack" << endl;
00769         s << "exch sub { end } repeat" << endl;
00770         s << "restore" << endl;
00771         s << "%%EOF" << endl;
00772     }
00773 
00774     /** Quadratic triangle shape function for use in write_eps_triangle */
00775     static double tN0( double r, double s )
00776     {
00777         double t = 1 - r - s;
00778         return t * ( 2 * t - 1 );
00779     }
00780     static double tN1( double r, double )
00781     {
00782         return r * ( 2 * r - 1 );
00783     }
00784     static double tN2( double, double s )
00785     {
00786         return s * ( 2 * s - 1 );
00787     }
00788     static double tN3( double r, double s )
00789     {
00790         double t = 1 - r - s;
00791         return 4 * r * t;
00792     }
00793     static double tN4( double r, double s )
00794     {
00795         return 4 * r * s;
00796     }
00797     static double tN5( double r, double s )
00798     {
00799         double t = 1 - r - s;
00800         return 4 * s * t;
00801     }
00802     /** Quadratic triangle shape function for use in write_eps_triangle */
00803     static Vector3D quad_tri_pt( double r, double s, const Vector3D* coords )
00804     {
00805         Vector3D result = tN0( r, s ) * coords[0];
00806         result += tN1( r, s ) * coords[1];
00807         result += tN2( r, s ) * coords[2];
00808         result += tN3( r, s ) * coords[3];
00809         result += tN4( r, s ) * coords[4];
00810         result += tN5( r, s ) * coords[5];
00811         return result;
00812     }
00813 
00814     void write_eps_triangle( Mesh* mesh, Mesh::ElementHandle elem, const char* filename, bool draw_iso_lines,
00815                              bool draw_nodes, MsqError& err, int width, int height )
00816     {
00817         // Get triangle vertices
00818         MsqVertex coords[6];
00819         EntityTopology type;
00820         mesh->elements_get_topologies( &elem, &type, 1, err );MSQ_ERRRTN( err );
00821         if( type != TRIANGLE )
00822         {
00823             MSQ_SETERR( err )( "Invalid element type", MsqError::UNSUPPORTED_ELEMENT );
00824             return;
00825         }
00826         std::vector< Mesh::VertexHandle > verts;
00827         std::vector< size_t > junk;
00828         mesh->elements_get_attached_vertices( &elem, 1, verts, junk, err );MSQ_ERRRTN( err );
00829         if( verts.size() != 3 && verts.size() != 6 )
00830         {
00831             MSQ_SETERR( err )( "Invalid element type", MsqError::UNSUPPORTED_ELEMENT );
00832             return;
00833         }
00834         mesh->vertices_get_coordinates( arrptr( verts ), coords, verts.size(), err );MSQ_ERRRTN( err );
00835 
00836         Vector3D coords2[6];
00837         std::copy( coords, coords + verts.size(), coords2 );
00838 
00839         std::vector< bool > fixed( verts.size(), false );
00840         if( draw_nodes )
00841         {
00842             mesh->vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err );MSQ_ERRRTN( err );
00843         }
00844         write_eps_triangle( coords2, verts.size(), filename, draw_iso_lines, draw_nodes, err, fixed, width, height );
00845     }
00846 
00847     void write_eps_triangle( const Vector3D* coords, size_t num_vtx, const char* filename, bool draw_iso_lines,
00848                              bool draw_nodes, MsqError& err, const std::vector< bool >& fixed, int width, int height )
00849     {
00850         const int PT_RAD        = 3;           // radius of circles for drawing nodes, in points
00851         const int PAD           = PT_RAD + 2;  // margin in points
00852         const double EDGE_GRAY  = 0.0;         // color for triangle edges, 0.0 => black
00853         const double ISO_GRAY   = 0.7;         // color for parameter iso-lines
00854         const double NODE_GRAY  = 0.0;         // color for node circle
00855         const double FIXED_GRAY = 1.0;         // color to fill fixed nodes with, 1.0 => white
00856         const double FREE_GRAY  = 0.0;         // color to fill free nodes with
00857 
00858         Projection proj( X, Y );
00859         Transform2D transf( coords, num_vtx, proj, width, height );
00860 
00861         // Open the file
00862         ofstream str( filename );
00863         if( !str )
00864         {
00865             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00866             return;
00867         }
00868 
00869         // Write header
00870         str << "%!PS-Adobe-2.0 EPSF-2.0" << endl;
00871         str << "%%Creator: Mesquite" << endl;
00872         str << "%%Title: Mesquite " << endl;
00873         str << "%%DocumentData: Clean7Bit" << endl;
00874         str << "%%Origin: 0 0" << endl;
00875         str << "%%BoundingBox: " << -PAD << ' ' << -PAD << ' ' << transf.max_horizontal() + PAD << ' '
00876             << transf.max_vertical() + PAD << endl;
00877         str << "%%Pages: 1" << endl;
00878 
00879         str << "%%BeginProlog" << endl;
00880         str << "save" << endl;
00881         str << "countdictstack" << endl;
00882         str << "mark" << endl;
00883         str << "newpath" << endl;
00884         str << "/showpage {} def" << endl;
00885         str << "/setpagedevice {pop} def" << endl;
00886         str << "%%EndProlog" << endl;
00887 
00888         str << "%%Page: 1 1" << endl;
00889         str << "1 setlinewidth" << endl;
00890         str << EDGE_GRAY << " setgray" << endl;
00891 
00892         const double h = 0.5, t = 1.0 / 3.0, w = 2.0 / 3.0, s = 1. / 6, f = 5. / 6;
00893         const int NUM_ISO                      = 15;
00894         const double iso_params[NUM_ISO][2][2] = {
00895             { { h, 0 }, { h, h } },  // r = 1/2
00896             { { t, 0 }, { t, w } },  // r = 1/3
00897             { { w, 0 }, { w, t } },  // r = 2/3
00898             { { s, 0 }, { s, f } },  // r = 1/6
00899             { { f, 0 }, { f, s } },  // r = 5/6
00900             { { 0, h }, { h, h } },  // s = 1/2
00901             { { 0, t }, { w, t } },  // s = 1/3
00902             { { 0, w }, { t, w } },  // s = 2/3
00903             { { 0, s }, { f, s } },  // s = 1/6
00904             { { 0, f }, { s, f } },  // s = 5/6
00905             { { 0, h }, { h, 0 } },  // t = 1 - r - s = 1/2
00906             { { 0, w }, { w, 0 } },  // t = 1 - r - s = 1/3
00907             { { 0, t }, { t, 0 } },  // t = 1 - r - s = 2/3
00908             { { 0, f }, { f, 0 } },  // t = 1 - r - s = 1/6
00909             { { 0, s }, { s, 0 } }   // t = 1 - r - s = 5/6
00910         };
00911 
00912         if( num_vtx == 3 )
00913         {
00914             int x[3], y[3];
00915             for( size_t i = 0; i < 3; ++i )
00916                 transf.transform( coords[i], x[i], y[i] );
00917 
00918             str << "newpath" << endl;
00919             str << x[0] << ' ' << y[0] << " moveto" << endl;
00920             str << x[1] << ' ' << y[1] << " lineto" << endl;
00921             str << x[2] << ' ' << y[2] << " lineto" << endl;
00922             str << x[0] << ' ' << y[0] << " lineto" << endl;
00923             str << "stroke" << endl;
00924 
00925             if( draw_iso_lines )
00926             {
00927                 str << ISO_GRAY << " setgray" << endl;
00928                 str << "newpath" << endl;
00929                 for( int i = 0; i < NUM_ISO; ++i )
00930                 {
00931                     double R[2]   = { iso_params[i][0][0], iso_params[i][1][0] };
00932                     double S[2]   = { iso_params[i][0][1], iso_params[i][1][1] };
00933                     double T[2]   = { 1 - R[0] - S[0], 1 - R[1] - S[1] };
00934                     Vector3D p[2] = { T[0] * coords[0] + R[0] * coords[1] + S[0] * coords[2],
00935                                       T[1] * coords[0] + R[1] * coords[1] + S[1] * coords[2] };
00936                     transf.transform( p[0], x[0], y[0] );
00937                     transf.transform( p[1], x[1], y[1] );
00938                     str << x[0] << ' ' << y[0] << " moveto" << endl;
00939                     str << x[1] << ' ' << y[1] << " lineto" << endl;
00940                 }
00941                 str << "    stroke" << endl;
00942             }
00943         }
00944         else if( num_vtx == 6 )
00945         {
00946             str << "newpath" << endl;
00947             write_eps_quadratic_edge( str, transf, coords[0], coords[3], coords[1] );
00948             write_eps_quadratic_edge( str, transf, coords[1], coords[4], coords[2] );
00949             write_eps_quadratic_edge( str, transf, coords[2], coords[5], coords[0] );
00950             str << "stroke" << endl;
00951 
00952             if( draw_iso_lines )
00953             {
00954                 str << ISO_GRAY << " setgray" << endl;
00955                 str << "newpath" << endl;
00956                 for( int i = 0; i < NUM_ISO; ++i )
00957                 {
00958                     double R[3]   = { iso_params[i][0][0], 0, iso_params[i][1][0] };
00959                     double S[3]   = { iso_params[i][0][1], 0, iso_params[i][1][1] };
00960                     R[1]          = 0.5 * ( R[0] + R[2] );
00961                     S[1]          = 0.5 * ( S[0] + S[2] );
00962                     Vector3D p[3] = { quad_tri_pt( R[0], S[0], coords ), quad_tri_pt( R[1], S[1], coords ),
00963                                       quad_tri_pt( R[2], S[2], coords ) };
00964                     write_eps_quadratic_edge( str, transf, p[0], p[1], p[2] );
00965                 }
00966                 str << "    stroke" << endl;
00967             }
00968         }
00969 
00970         if( draw_nodes )
00971         {
00972             for( size_t i = 0; i < num_vtx; ++i )
00973             {
00974                 int iw, ih;
00975                 // fill interior with either white or black depending
00976                 // on whether or not the vertex is fixed.
00977                 if( fixed[i] )
00978                     str << FIXED_GRAY << " setgray" << endl;
00979                 else
00980                     str << FREE_GRAY << " setgray" << endl;
00981                 transf.transform( coords[i], iw, ih );
00982                 str << w + PT_RAD << ' ' << ih << " moveto" << endl;
00983                 str << iw << ' ' << ih << ' ' << PT_RAD << " 0 360 arc" << endl;
00984                 str << "closepath fill" << endl;
00985                 str << NODE_GRAY << " setgray" << endl;
00986                 str << "newpath" << endl;
00987                 str << iw << ' ' << ih << ' ' << PT_RAD << " 0 360 arc" << endl;
00988                 str << "stroke" << endl;
00989             }
00990         }
00991 
00992         // Write footer
00993         str << "%%Trailer" << endl;
00994         str << "cleartomark" << endl;
00995         str << "countdictstack" << endl;
00996         str << "exch sub { end } repeat" << endl;
00997         str << "restore" << endl;
00998         str << "%%EOF" << endl;
00999     }
01000 
01001     void write_svg( Mesh* mesh, const char* filename, Projection proj, MsqError& err )
01002     {
01003         // Get a global patch
01004         PatchData pd;
01005         pd.set_mesh( mesh );
01006         pd.fill_global_patch( err );MSQ_ERRRTN( err );
01007 
01008         Transform2D transf( &pd, proj, 400, 400, true );
01009 
01010         // Open the file
01011         ofstream file( filename );
01012         if( !file )
01013         {
01014             MSQ_SETERR( err )( MsqError::FILE_ACCESS );
01015             return;
01016         }
01017 
01018         // Write header
01019         file << "<?xml version=\"1.0\" standalone=\"no\"?>" << endl;
01020         file << "<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\" "
01021              << "\"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">" << endl;
01022         file << endl;
01023         file << "<svg width=\"" << transf.max_horizontal() << "\" height=\"" << transf.max_vertical()
01024              << "\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" << endl;
01025 
01026         // Write mesh edges
01027         EdgeIterator iter( &pd, err );MSQ_ERRRTN( err );
01028         while( !iter.is_at_end() )
01029         {
01030             int s_w, s_h, e_w, e_h;
01031             transf.transform( iter.start(), s_w, s_h );
01032             transf.transform( iter.end(), e_w, e_h );
01033 
01034             file << "<line "
01035                  << "x1=\"" << s_w << "\" "
01036                  << "y1=\"" << s_h << "\" "
01037                  << "x2=\"" << e_w << "\" "
01038                  << "y2=\"" << e_h << "\" "
01039                  << " style=\"stroke:rgb(99,99,99);stroke-width:2\""
01040                  << "/>" << endl;
01041 
01042             iter.step( err );MSQ_ERRRTN( err );
01043         }
01044 
01045         // Write footer
01046         file << "</svg>" << endl;
01047     }
01048 
01049 }  // namespace MeshWriter
01050 
01051 }  // namespace MBMesquite
01052 
01053 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines