meshio.C

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include <tumble-conf.h>
#endif /* HAVE_CONFIG_H */

#include <fstream>
#include <cstdio>

#include "meshio.h"

#include "beziermesh.h"
#include "boundarymesh.h"
#include "hash_map.h"
#include "simulation.h"
#include "spline.h"

using namespace std;
using namespace hashers;

bool MeshInput::read(const char *filename) {
    char node_filename[ MAX_FILENAME ];
    char ele_filename[  MAX_FILENAME ];
    char edge_filename[ MAX_FILENAME ];
    char bdry_filename[ MAX_FILENAME ];

    bzero( node_filename, MAX_FILENAME );
    bzero( ele_filename,  MAX_FILENAME );
    bzero( edge_filename, MAX_FILENAME );
    bzero( bdry_filename, MAX_FILENAME );

    snprintf( node_filename, sizeof(node_filename), "%s.node", filename );
    snprintf( ele_filename,  sizeof(ele_filename),  "%s.ele",  filename );
    snprintf( edge_filename, sizeof(edge_filename), "%s.edge", filename );
    snprintf( bdry_filename, sizeof(bdry_filename), "%s.bdry", filename );

    node_file = fopen( node_filename, "r" );
    ele_file  = fopen( ele_filename, "r"  );
    edge_file = fopen( edge_filename, "r" );
    bdry_file = fopen( bdry_filename, "r" );

    if(!node_file) cout<<"MeshInput --- Error opening file: "<<node_filename<<endl;
    if(!ele_file)  cout<<"MeshInput --- Error opening file: "<< ele_filename<<endl;
    if(!edge_file) cout<<"MeshInput --- Error opening file: "<<edge_filename<<endl;
    if(!bdry_file) cout<<"MeshInput --- Error opening file: "<<bdry_filename<<endl;

    if( !node_file || !ele_file || !edge_file || !bdry_file ){
        if( fclose( node_file ) ) cout<<"Error closing node file!"<<endl;
        if( fclose( ele_file  ) ) cout<<"Error closing ele file!"<<endl;
        if( fclose( edge_file ) ) cout<<"Error closing edge file!"<<endl;
        if( fclose( bdry_file ) ) cout<<"Error closing bdry file!"<<endl;
        return false;
    }

    /* Step 1: read in node_file to the Node structure
     */
    vector<Node> nodes;
    char file_line[ MAX_FILE_LINE ];
    bzero( file_line, MAX_FILE_LINE );

    while ( fgets( file_line, MAX_FILE_LINE, node_file ) != NULL ) {
        nodes.push_back( Node( file_line ) );
    }
    //cout<<"Read "<<nodes.size()<<" Nodes"<<endl;
    /* Step 2: read BoundaryVertex's from bdry_file and add them to the
     * BoundaryMesh.  Only add CP's, we will add DP's later.
     */
    BoundaryVertexData bdry_vert_data;
    while ( read_boundary_vertex_data( &bdry_vert_data ) ) {
        bdry_verts[ bdry_vert_data.idx ] = bdry_mesh->add_boundary_vertex( nodes[ bdry_vert_data.node ].p );
    }
    //cout<<"Read "<<bdry_verts.size()<<" Boundary Vertexs"<<endl;
    /* Step 3: read in BoundaryEdge's from bdry_file.  Here we will only be adding the CP's.
     * The DP's will be added in later with the BezierEdge's.
     */
    BoundaryEdgeData bdry_edge_data;
    while ( read_boundary_edge_data( &bdry_edge_data ) ) {
        bdry_edges[ bdry_edge_data.idx ] = bdry_mesh->add_boundary_edge( bdry_verts[ bdry_edge_data.v0 ],
                                                                         bdry_verts[ bdry_edge_data.v1 ],
                                                                         bdry_edge_data.d,
                                                                         bdry_edge_data.k,
                                                                         (Movement)bdry_edge_data.fixed,
                                                                         bdry_edge_data.color,
                                                                         bdry_edge_data.restlength );
    }
    //cout<<"Read "<<bdry_edges.size()<<" Boundary Edges"<<endl;
    /* Step 4: Read in the BoundaryFace's and add them in */
    BoundaryFaceData bdry_face_data;
    bdry_faces[0] = NULL;

    while ( read_boundary_face_data( &bdry_face_data ) ) {
        bdry_faces[ bdry_face_data.idx ] = bdry_mesh->add_boundary_face( bdry_face_data.edges,
                                                                         bdry_face_data.min_angle,
                                                                         bdry_face_data.color );
    }
    //cout<<"Read "<<bdry_faces.size()<<" Boundary Faces"<<endl;
    /* Step 5: Read in the BezierEdge's from edge_file.  For each edge we will first consider
     * its vertexs.  If these coorespond to nodes in the node list that have not yet been added
     * then we will add the cooresponding BezierVertex along with its proper CP and DP.
     *
     * However there is an exception to be made in the case of BezierEdge's on the boundary.
     * Assuming the coorespoing nodes have not been added, we must perform a reverse lookup on the
     * cooresponding BoundaryEdge, via the u0 and u1 values.  This will give us a pointer to the CP.
     * we need to use this to hang the DP as we skiped that before.  Also we will need to hang data
     * off the internal point, whose CP was created by the QBSpline constructor in Step 3.
     *
     * Finally when all these points have been added or updated, we may add in the Edge.
     * Additionally, we must be sure to add all boundary edges BEFORE all non-boundary edges, thus we
     * keep a vector of the non-boundary edges and deal with them later.
     */
    BezierEdgeData edge_data;
    Node *node;
    vector<BezierEdgeData> internal_edges;
    ControlPoint cp0, cp1, cp2;
    int cp0_is_vert, cp2_is_vert;
    BoundaryEdge *be;
    BezierEdge *edge;
    BezierVertex *vert;
    while ( read_bezier_edge_data( &edge_data ) ) {
        if ( edge_data.is_bdry ) {

            be = bdry_edges[ edge_data.bdry ];
            if ( ! be->get_spline()->get_edge_cps( edge_data.u0, edge_data.u1, cp0, cp1, cp2, cp0_is_vert, cp2_is_vert ) ) {
                printf( "Simulation::read_from_file  error, could not get edge cps from boundary mesh\n" );
                return false;
            }

            node = &nodes[ edge_data.v0 ];
            if ( !node->added ) {
                node->added = true;
                vert = bezier_mesh->add_bezier_vertex( cp0, node->d );
                bez_verts[ edge_data.v0 ] = vert;
                if ( cp0_is_vert == 0 )
                    vert->set_bdry( be->get_vertex(0) );
                else if ( cp0_is_vert == 1 )
                    vert->set_bdry( be->get_vertex(1) );
                else
                    vert->set_bdry( be, edge_data.u0 );
            }

            node = &nodes[ edge_data.v2 ];
            if ( !node->added ) {
                node->added = true;
                vert = bezier_mesh->add_bezier_vertex( cp2, node->d );
                bez_verts[ edge_data.v2 ] = vert;
                if ( cp2_is_vert == 0 )
                    vert->set_bdry( be->get_vertex(0) );
                else if ( cp2_is_vert == 1 )
                    vert->set_bdry( be->get_vertex(1) );
                else
                    vert->set_bdry( be, edge_data.u1 );
            }

            node = &nodes[ edge_data.v1 ];
            node->added = true;
            edge = bezier_mesh->add_bezier_edge( cp1, node->d, bez_verts[ edge_data.v0 ], bez_verts[ edge_data.v2 ] );
            bez_edges[ edge_data.idx ] = edge;
            edge->set_bdry( be, edge_data.u0, edge_data.u1 );
        } else {
            internal_edges.push_back( edge_data );
        }


    }
    //cout<<"Read "<<bez_verts.size()<<" Bezier Vertexs"<<endl;

    /* Step 5b: Now add all internal edges */
    for ( vector<BezierEdgeData>::iterator i = internal_edges.begin(); i != internal_edges.end(); i++ ) {

        node = &nodes[ i->v0 ];
        if ( !node->added ) {
            bez_verts[ i->v0 ] = bezier_mesh->add_bezier_vertex( node->p, node->d );
            node->added = true;
        }
        node = &nodes[ i->v2 ];
        if ( !node->added ) {
            bez_verts[ i->v2 ] = bezier_mesh->add_bezier_vertex( node->p, node->d );
            node->added = true;
        }
        node = &nodes[ i->v1 ];
        bez_edges[ i->idx ] = bezier_mesh->add_bezier_edge( node->p, node->d, bez_verts[ i->v0 ], bez_verts[ i->v2 ] );
    }
    //cout<<"Read "<<bez_edges.size()<<" Bezier Edges"<<endl;

    /* Step 6: Read in the BezierTriangle's from ele_file.
     */
    BezierTriangleData tri_data;
    int num_bez_tris=0;
    while ( read_bezier_triangle_data( &tri_data ) ) {
        bezier_mesh->add_bezier_triangle( bdry_faces[ tri_data.bdry_face ], bez_edges[ tri_data.e0 ], bez_edges[ tri_data.e1 ], bez_edges[ tri_data.e2 ] );
        num_bez_tris++;
    }
    //cout<<"Read "<<num_bez_tris<<" Bezier Triangles"<<endl;

    if( fclose( node_file ) ) cout<<"Error closing node file!"<<endl;
    if( fclose( ele_file  ) ) cout<<"Error closing ele file!"<<endl;
    if( fclose( edge_file ) ) cout<<"Error closing edge file!"<<endl;
    if( fclose( bdry_file ) ) cout<<"Error closing bdry file!"<<endl;

    return true;
}

MeshInput::Node::Node( char *line ) {
    char * str;
    vector<double> data;
    added = FALSE;
    str = strtok( line, " \n" );
    num = atoi( str );
    str = strtok( NULL, " \n" );
    p.coords[ 0 ] = strtod( str, NULL );
    str = strtok( NULL, " \n" );
    p.coords[ 1 ] = strtod( str, NULL );
    str = strtok( NULL, " \n" );
    while ( str != NULL ) {
        data.push_back( strtod( str, NULL ) );
        str = strtok( NULL, " \n" );
    }
    d.update( data );
}

void MeshInput::Node::print() {
    cout << "Node: num: " << num << " added: " << added << " p: " << p << " d: " << d << endl;
}

int MeshInput::read_boundary_vertex_data( BoundaryVertexData *data ) {
    char str[ MAX_FILE_LINE ];
    char *temp;
    if ( ( fgets( str, MAX_FILE_LINE, bdry_file ) == NULL ) || !strncmp( str, "---", 3 ) ) {
        return FALSE;
    }
    temp = strtok( str, " \n" );
    data->idx = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->node = atoi( temp );
    return TRUE;
}

int MeshInput::read_boundary_edge_data( BoundaryEdgeData *data ) {
    const int max_size = 10*MAX_FILE_LINE;
    char str[ max_size ];
    char *temp;
    int closed;
    int num;
    float f;
    Point2D d;
    char *ret;

    ret = fgets( str, max_size, bdry_file );
    if( (ret == NULL) || (strlen(str) >= max_size-2) ){
        cout<<"MeshInput::read_boundary_edge_data -- error reading first line or line too long"<<endl;
        return FALSE;
    }
    if ( !strncmp( str, "===", 3 ) ) {
        /* done with boundary edges */
        return FALSE;
    }
    num = sscanf( str, "%u %u %u %i %i %i %f", &data->idx, &data->v0, &data->v1, &closed, &data->fixed , &data->color, &f );
    if(num == 6)        data->restlength=0.0;
    else if(num == 7)   data->restlength=(double)f;
    else {
    printf("MeshInput::read_boundary_edge_data --- error reading spline data\n");
    exit(0);
    }

    closed ? data->closed = true : data->closed = false;

    data->d.clear();
    data->k.clear();
    /* Read in next line with d values */
    ret=fgets( str, max_size, bdry_file );
    if (( ret == NULL ) || (strlen(str)>= max_size-2)){
        cout<<"MeshInput::read_boundary_edge_data -- error reading second line or line too long"<<endl;
        return FALSE;
    }

    temp = strtok( str, " \n" );
    while ( temp != NULL ) {
        d.coords[ 0 ] = strtod( temp, NULL );
        temp = strtok( NULL, " \n" );
        d.coords[ 1 ] = strtod( temp, NULL );
        data->d.push_back( d );
        temp = strtok( NULL, " \n" );
    }

    /* Read in last line with k values */
    ret=fgets( str, max_size, bdry_file );
    if (( ret == NULL ) || (strlen(str) >= max_size-2)){
        cout<<"MeshInput::read_boundary_edge_data -- error reading third line or line too long"<<endl;
        return FALSE;
    }

    temp = strtok( str, " \n" );
    while ( temp != NULL ) {
        data->k.push_back( strtod( temp, NULL ) );
        temp = strtok( NULL, " \n" );
    }
    return TRUE;
}

int MeshInput::read_boundary_face_data( BoundaryFaceData *data ) {
    char str[ MAX_FILE_LINE ];
    char *temp;

    if ( fgets( str, MAX_FILE_LINE, bdry_file ) == NULL ) {
        return FALSE;
    }
    data->edges.clear();
    temp = strtok( str, " \n" );
    data->idx = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->min_angle = strtod( temp, NULL );
    temp = strtok( NULL, " \n" );
    data->color = atoi( temp );
    while ( ( temp = strtok( NULL, " \n" ) ) != NULL ) {
        data->edges.push_back( bdry_edges[ atoi( temp ) ] );
    }
    return TRUE;
}

int MeshInput::read_bezier_edge_data( BezierEdgeData *data ) {
    char str[ MAX_FILE_LINE ];
    bzero( str, MAX_FILE_LINE );
    char *temp;
    if ( fgets( str, MAX_FILE_LINE, edge_file ) == NULL ) {
        return FALSE;
    }
    temp = strtok( str, " \n" );
    data->idx = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->v0 = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->v1 = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->v2 = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->bdry = atoi( temp );
    ( data->bdry < 0 ) ? data->is_bdry = false : data->is_bdry = true;

    if ( data->is_bdry ) {
        temp = strtok( NULL, " \n" );
        data->u0 = strtod( temp, NULL );
        temp = strtok( NULL, " \n" );
        data->u1 = strtod( temp, NULL );
    } else {
        data->u0 = 0.0;
        data->u1 = 0.0;
    }
    return TRUE;
}

int MeshInput::read_bezier_triangle_data( BezierTriangleData *data ) {
    char str[ MAX_FILE_LINE ];
    char *temp;
    int inverted;
    if ( fgets( str, MAX_FILE_LINE, ele_file ) == NULL ) {
        return FALSE;
    }
    temp = strtok( str, " \n" );
    data->idx = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->e0 = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->e1 = atoi( temp );
    temp = strtok( NULL, " \n" );
    data->e2 = atoi( temp );
    temp = strtok( NULL, " \n" );
    inverted = atoi( temp );
    inverted ? data->inverted = true : data->inverted = false;
    temp = strtok( NULL, " \n" );
    data->bdry_face = atoi( temp );
    return TRUE;
}

bool MeshOutput::write(const char *filename) {
    char node_filename[ MAX_FILENAME ];
    char ele_filename[ MAX_FILENAME ];
    char edge_filename[ MAX_FILENAME ];
    char bdry_filename[ MAX_FILENAME ];

    snprintf( node_filename, sizeof(node_filename), "%s.node", filename );
    snprintf( ele_filename,  sizeof(ele_filename),  "%s.ele",  filename );
    snprintf( edge_filename, sizeof(edge_filename), "%s.edge", filename );
    snprintf( bdry_filename, sizeof(bdry_filename), "%s.bdry", filename );

    node_file = fopen( node_filename, "w" );
    ele_file  = fopen( ele_filename, "w"  );
    edge_file = fopen( edge_filename, "w" );
    bdry_file = fopen( bdry_filename, "w" );

    if(!node_file) cout<<"MeshInput --- Error opening file: "<<node_filename<<endl;
    if(!ele_file)  cout<<"MeshInput --- Error opening file: "<< ele_filename<<endl;
    if(!edge_file) cout<<"MeshInput --- Error opening file: "<<edge_filename<<endl;
    if(!bdry_file) cout<<"MeshInput --- Error opening file: "<<bdry_filename<<endl;

    if(!node_file || !ele_file || !edge_file || !bdry_file){
        if( fclose( node_file ) ) cout<<"Error closing node file!"<<endl;
        if( fclose( ele_file  ) ) cout<<"Error closing ele file!"<<endl;
        if( fclose( edge_file ) ) cout<<"Error closing edge file!"<<endl;
        if( fclose( bdry_file ) ) cout<<"Error closing bdry file!"<<endl;
        return false;
    }


    /* These integers are used to assign an enumeration to control points and cells via
     *the hash_maps */
    unsigned node_index = 0;
    unsigned bez_edge_index = 0;
    unsigned bdry_vert_index = 0;
    unsigned bdry_edge_index = 0;
    unsigned bdry_face_index = 0;

    /* Step 1: dump all control points to node file.  This could be done by iterating
     * through the list of ControlPoints, however this will not let us determine
     * which data points are associated with each control point, therefore we
     * will iterate through the Bezier Vertexs.  I am assuming that BezierVertexs,
     * ControlPoints, and DataPoints are all in one-to-one coorespondence.
     */
    BezierMesh::Vertex_Hash_T::iterator viter;
    for(viter = bezier_mesh->get_vertices_begin();
        viter != bezier_mesh->get_vertices_end(); ++viter) {
        BezierVertex *v = *viter;
        node_hash[v->get_control_point()] = node_index++;
        node_to_file( v->get_control_point(), v->get_data_point() );
    }

    /* Step 2: dump all boundary vertexs to file*/
    BoundaryMesh::Vertex_Hash_T::iterator bviter;
    for(bviter = bdry_mesh->get_vertices_begin();
        bviter != bdry_mesh->get_vertices_end(); ++bviter ) {
        BoundaryVertex *bv = *bviter;
        bdry_vert_hash[ bv ] = bdry_vert_index++;
        bdry_vert_to_file( bv );
    }
    fprintf( bdry_file, "--- %u\n", bdry_vert_index );

    /* Step 3: dump all boundary edges to file*/
    BoundaryMesh::Edge_Hash_T::iterator beiter;
    for( beiter = bdry_mesh->get_edges_begin();
        beiter != bdry_mesh->get_edges_end(); ++beiter ) {
        BoundaryEdge *be = *beiter;
        bdry_edge_hash[ be ] = bdry_edge_index++;
        bdry_edge_to_file( be );
    }
    fprintf( bdry_file, "=== %u\n", bdry_edge_index );

    /* Step 4: dump all boundary faces to file*/
    BoundaryMesh::Face_Hash_T::iterator bfiter;
    for( bfiter = bdry_mesh->get_faces_begin();
        bfiter != bdry_mesh->get_faces_end(); ++bfiter ) {
        BoundaryFace *bf = *bfiter;
        bdry_face_hash[ bf ] = bdry_face_index++;
        bdry_face_to_file( bf );
    }

    /* Step 5: dump all bezier edges to file, but we also need to add thier
     *internal control point to the node file
     */
    BezierMesh::Edge_Hash_T::iterator eiter;
    for(eiter = bezier_mesh->get_edges_begin();
        eiter != bezier_mesh->get_edges_end(); ++eiter) {
        BezierEdge *e = *eiter;
        bez_edge_hash[ e ] = bez_edge_index++;
        node_hash[ e->get_control_point(1) ] = node_index++;
        node_to_file( e->get_control_point(1), e->get_data_point(1) );
        edge_to_file( e );
    }

    /* Step 6: dump all bezier triangles to file.  Here we do not need a hash because
     * these are the last to be added.  Instead the ele_to_file() function holds a static
     * integer index for the triangles to number them
     */
    BezierMesh::Face_Hash_T::iterator fiter;
    for(fiter = bezier_mesh->get_faces_begin();
        fiter != bezier_mesh->get_faces_end(); ++fiter ) {
        BezierTriangle *f = *fiter;
        ele_to_file( f, bezier_mesh);
    }

    if( fclose( node_file ) ) cout<<"Error closing node file!"<<endl;
    if( fclose( ele_file  ) ) cout<<"Error closing ele file!"<<endl;
    if( fclose( edge_file ) ) cout<<"Error closing edge file!"<<endl;
    if( fclose( bdry_file ) ) cout<<"Error closing bdry file!"<<endl;
    return true;
}


int MeshOutput::node_to_file( ControlPoint cp, DataPoint dp ) {
    fprintf( node_file, "%u %.12f %.12f", node_hash[ cp ], cp->coords[ 0 ], cp->coords[ 1 ] );
    for ( int i = 0; i < dp->length; i++ ) {
        fprintf( node_file, " %.12f", dp->fields[ i ] );
    }
    fprintf( node_file, "\n" );
    return TRUE;
}

int MeshOutput::edge_to_file( BezierEdge* edge ) {
    if ( ! edge->is_boundary() ) {
        fprintf( edge_file, "%u %u %u %u -1\n",
                 bez_edge_hash[ edge ],
                 node_hash[ edge->get_control_point(0) ],
                 node_hash[ edge->get_control_point(1) ],
                 node_hash[ edge->get_control_point(2) ] );
    } else {
        fprintf( edge_file, "%u %u %u %u %u %.12f %.12f\n",
                 bez_edge_hash[edge],
                 node_hash[ edge->get_control_point(0) ],
                 node_hash[ edge->get_control_point(1) ],
                 node_hash[ edge->get_control_point(2) ],
                 bdry_edge_hash[ edge->get_bdry_edge() ],
                 edge->get_u(0), edge->get_u(1) );
    }
    return TRUE;
}

int MeshOutput::ele_to_file( BezierTriangle* tri, BezierMesh *bezier_mesh ) {
    static unsigned tri_number = 0;
    fprintf( ele_file, "%u %u %u %u %u %u\n", tri_number++,
             bez_edge_hash[ tri->get_edge(0) ],
             bez_edge_hash[ tri->get_edge(1) ],
             bez_edge_hash[ tri->get_edge(2) ],
             tri->is_inverted( ) ,
             bdry_face_hash[ tri->get_bdry_face_or_null() ] );
    return TRUE;
}

int MeshOutput::bdry_vert_to_file( BoundaryVertex *vert ) {
    fprintf( bdry_file, "%u %u\n",
        bdry_vert_hash[ vert ], node_hash[ vert->get_control_point() ] );
    return TRUE;
}

int MeshOutput::bdry_edge_to_file( BoundaryEdge *edge ) {
    fprintf( bdry_file, "%u %u %u %u %u %u %.12f\n", bdry_edge_hash[ edge ],
             bdry_vert_hash[edge->get_vertex(0)],
             bdry_vert_hash[edge->get_vertex(1)],
             edge->get_spline()->closed, edge->get_fixed(), edge->get_color(),
             edge->get_spline()->restlength )
    ;
    int num_deboor = edge->get_spline()->get_num_deboor();
    ControlPoint d;
    unsigned i;
    for ( i=0; i < ( unsigned ) num_deboor; i++ ) {
        d = edge->get_spline()->get_deboor( i );
        fprintf( bdry_file, " %.12f %.12f", d->coords[ 0 ], d->coords[ 1 ] );
    }
    fprintf( bdry_file, "\n" );
    for ( i=0; i < edge->get_spline()->k.size(); i++ ) {
        fprintf( bdry_file, " %.12f", edge->get_spline()->k[ i ] );
    }
    fprintf( bdry_file, "\n" );
    return TRUE;
}

int MeshOutput::bdry_face_to_file( BoundaryFace *f ) {
    fprintf(bdry_file, "%u %.12f %u", bdry_face_hash[f],
            f->get_min_angle(), f->get_color() );
    for(BoundaryFace::edge_iterator it = f->begin_edges();
        it != f->end_edges(); ++it) {
      fprintf(bdry_file, " %u", bdry_edge_hash[*it]);
    }
    fprintf( bdry_file, "\n" );
    return TRUE;
}

static const char *official_header="Tumble Binary Geometry File";
static const int official_header_len = 27;
static const unsigned current_format_version=1;

MeshBinaryInput::MeshBinaryInput(Simulation* _sim)
{
    sim = _sim;
    bezier_mesh = sim->get_bezier_mesh();
    bdry_mesh = sim->get_boundary_mesh();
}

bool MeshBinaryInput::read(const char *filename)
{
    if(!filename) return false;
    int len = strlen(filename);
    if( strncmp(filename+(len-4), ".geo", 4) ){
        cout<<"MeshBinaryInput -- filename:[ \""<<filename<<"\" ]  does not end in extension \".geo\""<<endl;
        return false;
    }

    ifstream in( filename, ios::binary | ios::in );
    if(!in){
        cout<<"MeshBinaryInput -- Cannot open file for input: "<<filename<<endl;
        return false;
    }

    struct Node {
        Point2D p;
        LinearData d;
    };

    hash_map<unsigned, BoundaryVertex*> bdry_verts;
    hash_map<unsigned, BoundaryEdge*> bdry_edges;
    hash_map<unsigned, BoundaryFace*> bdry_faces;
    hash_map<unsigned, BezierVertex*> bez_verts;
    hash_map<unsigned, BezierEdge*> bez_edges;

    unsigned num_nodes;
    unsigned num_bdry_verts;
    unsigned num_bdry_edges;
    unsigned num_bdry_faces;
    unsigned num_bez_verts;
    unsigned num_bez_edges;
    unsigned num_bez_tris;

    unsigned i,j;
    unsigned node_num;
    unsigned v0_idx;
    unsigned v1_idx;
    unsigned e0_idx, e1_idx, e2_idx;
    unsigned e_num;
    unsigned f_num;
    unsigned fixed;
    unsigned color;

    //read in and check header
    char header[official_header_len+1];
    bzero(header,official_header_len+1);
    in.read(header, official_header_len);
    if(strncmp(official_header, header, official_header_len)){
        cout<<"MeshBinaryInput --- error recieved bad header: "<<header<<endl;
        return false;
    }

    //read in and check version number
    unsigned vers_num;
    in.read((char*) & vers_num, sizeof(unsigned));
    if( vers_num != current_format_version ){
        cout<<"MeshBinaryInput --- Attempting to read a binary file of incorrect version."<<endl
            <<"Current version: "<<current_format_version<<" File version: "<<vers_num<<endl;
        return false;
    }

    //read in nodes
    //recall indexing starts at 1, so reflect that in node list
    in.read((char*) &num_nodes, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_nodes<<" nodes"<<endl;
    Node *nodes = new Node[num_nodes+1];
    double *data = NULL;
    unsigned this_data_size;
    unsigned data_size = 0;
    for(i=1; i<num_nodes+1; i++){
        in.read( (char*) &node_num, sizeof(unsigned));
        in.read( (char*) &nodes[i].p, sizeof(Point2D));
        in.read( (char*) &this_data_size, sizeof(unsigned));
        assert(node_num == i);
        //assert(this_data_size < 100);
        if(data==NULL) {
            data_size=this_data_size;
            data = new double[data_size];
        }
        assert(data_size == this_data_size);
        for(j=0; j<data_size; j++){
            in.read((char*) &data[j], sizeof(double));
        }
        nodes[i].d = LinearData(data, data_size);
    }

    //read in boundary vertexs
    in.read((char*) &num_bdry_verts, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_bdry_verts<<" BoundaryVertexs"<<endl;
    unsigned bdry_vert_idx;
    for(i=0; i<num_bdry_verts; i++){
        in.read((char*) &bdry_vert_idx, sizeof(unsigned));
        in.read((char*) &node_num, sizeof(unsigned));
        in.read((char*) &fixed, sizeof(unsigned));
        assert(bdry_vert_idx == i+1);
        assert(node_num != 0 && node_num <= num_nodes);

        //create new BoundaryVertex
        //fixed isn't used as we can figure this out later from the boundary edges
        bdry_verts[bdry_vert_idx] = bdry_mesh->add_boundary_vertex( nodes[node_num].p );
    }

    //read in boundary edges
    in.read((char*) &num_bdry_edges, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_bdry_edges<<" BoundaryEdges"<<endl;
    unsigned bdry_edge_idx;
    unsigned num_internal_deboor;
    unsigned num_knots;
    double restlength;
    vector<Point2D> deboor;
    vector<double> knots;
    double kval;
    for(i=0; i<num_bdry_edges; i++){
        in.read((char*) &bdry_edge_idx, sizeof(unsigned));
        in.read((char*) &v0_idx, sizeof(unsigned));
        in.read((char*) &v1_idx, sizeof(unsigned));
        in.read((char*) &fixed, sizeof(unsigned));
        in.read((char*) &color, sizeof(unsigned));
        in.read((char*) &restlength, sizeof(double));

        assert(bdry_edge_idx == i+1);
        assert(v0_idx != 0 && v0_idx <= num_bdry_verts);
        assert(v1_idx != 0 && v1_idx <= num_bdry_verts);
        assert(restlength >= 0.0 );

        //read in deboor points
        //must add endpoints as well
        in.read((char*) &num_internal_deboor, sizeof(unsigned));
        deboor.clear();
        deboor.push_back(bdry_verts[v0_idx]->get_cp());
        for(j=0; j<num_internal_deboor; j++){
            in.read((char*) &node_num, sizeof(unsigned));
            deboor.push_back(nodes[node_num].p);
        }
        deboor.push_back(bdry_verts[v1_idx]->get_cp());

        //read in knot values
        knots.clear();
        in.read((char*) &num_knots, sizeof(unsigned));
        for(j=0; j<num_knots;j++){
            in.read((char*) &kval, sizeof(double));
            knots.push_back(kval);
        }
        assert(num_knots-1 == num_internal_deboor);
        assert(num_knots == deboor.size()-1);

        //create new BoundaryEdge
        bdry_edges[bdry_edge_idx] = bdry_mesh->add_boundary_edge( bdry_verts[v0_idx], bdry_verts[v1_idx],
                                                                  deboor, knots, (Movement)fixed, color, restlength);
    }

    // read in boundary faces
    in.read((char*) &num_bdry_faces, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_bdry_faces<<" BoundaryFaces"<<endl;
    unsigned bdry_face_idx;
    double min_angle;
    unsigned num_edges;
    vector<BoundaryEdge*> face_edges;
    for(i=0; i<num_bdry_faces; i++){
        in.read((char*) &bdry_face_idx, sizeof(unsigned));
        in.read((char*) &min_angle, sizeof(double));
        in.read((char*) &color, sizeof(unsigned));
        in.read((char*) &num_edges, sizeof(unsigned));
        assert(bdry_face_idx == i+1);
        assert(min_angle >= 0.0);
        assert(num_edges >= 1);

        face_edges.clear();
        for(j=0; j < num_edges; j++){
            in.read((char*) &e_num, sizeof(unsigned));
            face_edges.push_back( bdry_edges[e_num] );
        }

        bdry_faces[bdry_face_idx] = bdry_mesh->add_boundary_face( face_edges, min_angle, color);
    }

    //read in bezier vertexs
    in.read((char*) &num_bez_verts, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_bez_verts<<" BezierVertexs"<<endl;
    unsigned bez_vert_idx;
    unsigned bdry_type;
    unsigned bdry_idx;
    unsigned bdry_knot_idx;
    QBSpline *spline;
    BezierVertex *v;
    BoundaryVertex *bv;
    BoundaryEdge *be;
    for(i=0; i<num_bez_verts; i++){
        in.read((char*) &bez_vert_idx,  sizeof(unsigned));
        in.read((char*) &node_num,      sizeof(unsigned));
        in.read((char*) &bdry_type,     sizeof(unsigned));
        in.read((char*) &bdry_idx,      sizeof(unsigned));
        in.read((char*) &bdry_knot_idx, sizeof(unsigned));

        assert(bez_vert_idx == i+1);
        assert(node_num != 0 && node_num <= num_nodes);
        assert(bdry_type <= 2);

        if(bdry_type == 0){
            //d2 - internal vertex - no bdry to worry about
            assert(bdry_idx == 0);
            assert(bdry_knot_idx == 0);
            bez_verts[bez_vert_idx] = bezier_mesh->add_bezier_vertex( nodes[node_num].p, nodes[node_num].d );
        } else if(bdry_type==1){
            //d0 vertex - need to get CP from vertex and set bdry
            assert(bdry_idx != 0 && bdry_idx <= num_bdry_verts);
            assert(bdry_knot_idx == 0);
            bv=bdry_verts[bdry_idx];
            v = bezier_mesh->add_bezier_vertex( bv->get_control_point(),
                nodes[node_num].d );
            bez_verts[bez_vert_idx] = v;
            v->set_bdry(bv);
        } else{
            //d1 vertex - need to get CP from bdry edge with given knot idx and set bdry
            assert(bdry_idx != 0 && bdry_idx <= num_bdry_edges);
            assert(bdry_knot_idx != 0);
            be = bdry_edges[bdry_idx];
            spline = be->get_spline();
            v = bezier_mesh->add_bezier_vertex( spline->get_control_point_at_idx(bdry_knot_idx), nodes[node_num].d );
            bez_verts[bez_vert_idx] = v;
            v->set_bdry(be,(int)bdry_knot_idx);
        }
    }

    //read in bezier edges
    in.read((char*) &num_bez_edges, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_bez_edges<<" BezierEdges"<<endl;
    unsigned bez_edge_idx;
    BezierEdge *e;
    unsigned bdry_knot_idx0, bdry_knot_idx1;
    for(i=0; i<num_bez_edges; i++){
        in.read((char*) &bez_edge_idx,   sizeof(unsigned));
        in.read((char*) &v0_idx,         sizeof(unsigned));
        in.read((char*) &v1_idx,         sizeof(unsigned));
        in.read((char*) &node_num,       sizeof(unsigned));
        in.read((char*) &bdry_idx,       sizeof(unsigned));
        in.read((char*) &bdry_knot_idx0, sizeof(unsigned));
        in.read((char*) &bdry_knot_idx1, sizeof(unsigned));

        assert(bez_edge_idx == i+1);
        assert(v0_idx != 0 && v0_idx <= num_bez_verts);
        assert(v1_idx != 0 && v0_idx <= num_bez_verts);
        assert(node_num != 0 && node_num <= num_nodes);

        if(bdry_idx == 0){
            //not a bdry edge
            assert(bdry_knot_idx0 == 0);
            assert(bdry_knot_idx1 == 0);
            bez_edges[bez_edge_idx] = bezier_mesh->add_bezier_edge( nodes[node_num].p, nodes[node_num].d, bez_verts[v0_idx], bez_verts[v1_idx] );
        } else {
            //is a bdry edges, get node cp from spline, set bdry
            assert(bdry_idx <= num_bdry_edges);
            be = bdry_edges[bdry_idx];
            spline = be->get_spline();
            assert( abs((int)bdry_knot_idx0 - (int)bdry_knot_idx1) == 1 );
            e = bezier_mesh->add_bezier_edge( spline->get_edge_center_cp(bdry_knot_idx0, bdry_knot_idx1),
                                              nodes[node_num].d, bez_verts[v0_idx], bez_verts[v1_idx] );
            bez_edges[bez_edge_idx] = e;
            e->set_bdry(be, spline->get_u_at_idx(bdry_knot_idx0), spline->get_u_at_idx(bdry_knot_idx1));
        }
    }

    //read in bezier faces
    in.read((char*) &num_bez_tris, sizeof(unsigned));
    //cout<<"Attempting to read in "<<num_bez_tris<<" BezierTriangles"<<endl;
    for(i=0; i<num_bez_tris; i++){
        in.read((char*) &e0_idx, sizeof(unsigned));
        in.read((char*) &e1_idx, sizeof(unsigned));
        in.read((char*) &e2_idx, sizeof(unsigned));
    bool ihf; // Inverse-Handed-Face Flag  // tp517: 4/16/07
    in.read((char*) &ihf, sizeof(bool));
        in.read((char*) &f_num, sizeof(unsigned));

        assert(e0_idx != 0 && e0_idx <= num_bez_edges);
        assert(e1_idx != 0 && e1_idx <= num_bez_edges);
        assert(e2_idx != 0 && e2_idx <= num_bez_edges);
        /* f_num zero is the null face, so it's allowed */
        assert(f_num <= num_bdry_faces);

        bezier_mesh->add_bezier_triangle(bdry_faces[f_num], bez_edges[e0_idx], bez_edges[e1_idx], bez_edges[e2_idx],ihf); // read/write ihf: tp517 4/16/07
    }

    return true;
}

MeshBinaryOutput::MeshBinaryOutput(Simulation* _sim)
{
    sim = _sim;
    bezier_mesh = sim->get_bezier_mesh();
    bdry_mesh = sim->get_boundary_mesh();
}

MeshBinaryOutput::MeshBinaryOutput(BezierMesh *_bezier_mesh, BoundaryMesh *_bdry_mesh)
{
    sim = NULL;
    bezier_mesh = _bezier_mesh;
    bdry_mesh = _bdry_mesh;
}

bool MeshBinaryOutput::write(const char *filename)
{
    if(!filename) return false;
    int len = strlen(filename);
    if( strncmp(filename+(len-4), ".geo", 4) ){
        cout<<"MeshBinaryOutput -- filename:[ \""<<filename<<"\" ]  does not end in extension \".geo\""<<endl;
        return false;
    }

    ofstream out( filename, ios::binary | ios::out | ios::trunc );
    if(!out){
        cout<<"MeshBinaryOutput -- Cannot open file for output: "<<filename<<endl;
        return false;
    }

    /* write header and version number */
    out.write( official_header, official_header_len);
    out.write((const char *) &current_format_version, sizeof(unsigned));

    //iterators for cells
    BezierMesh::Vertex_Hash_T::const_iterator bez_v_it;
    BezierMesh::Edge_Hash_T::const_iterator bez_e_it;
    BezierMesh::Face_Hash_T::const_iterator bez_t_it;
    BoundaryMesh::Vertex_Hash_T::const_iterator bdry_v_it;
    BoundaryMesh::Edge_Hash_T::const_iterator bdry_e_it;
    BoundaryMesh::Face_Hash_T::const_iterator bdry_f_it;

    //mappings from cells to thier indexes
    hash_map<ControlPoint, unsigned, ControlPointHasher> node_idx;
    hash_map<BezierVertex*, unsigned> bezvert_idx;
    hash_map<BezierEdge*, unsigned> bezedge_idx;
    hash_map<BoundaryVertex*, unsigned> bdryvert_idx;
    hash_map<BoundaryEdge*, unsigned> bdryedge_idx;
    hash_map<BoundaryFace*, unsigned> bdryface_idx;

    //counters
    // we will start all indexing at 1, to provide 0 as a free variable with no meaning for special cases
    unsigned nodes_ctr = 1;
    unsigned bezverts_ctr = 1;
    unsigned bezedges_ctr = 1;
    unsigned beztris_ctr = 1;
    unsigned bdryverts_ctr = 1;
    unsigned bdryedges_ctr = 1;
    unsigned bdryfaces_ctr = 1;

    ControlPoint cp;
    DataPoint dp;
    unsigned data_len;
    unsigned node_num;
    unsigned v_num;
    unsigned e_num;
    unsigned f_num;
    unsigned num_edges;
    unsigned fixed;
    unsigned color;
    unsigned i;

    /* write data nodes.  These are all the data points and thier coordinates
     * there is one for every vertex and one for the middle cp of every edge
     */
    unsigned total_num_nodes = bezier_mesh->get_num_vertices()
                             + bezier_mesh->get_num_edges();
    out.write((char*) &total_num_nodes, sizeof(total_num_nodes));
    //bezier vert nodes
    for(bez_v_it = bezier_mesh->get_vertices_begin();
        bez_v_it != bezier_mesh->get_vertices_end(); ++bez_v_it){
        const BezierVertex *v = *bez_v_it;
        ControlPoint cp = v->get_control_point();
        node_idx[ cp ] = nodes_ctr;
        out.write((const char*) &nodes_ctr, sizeof(nodes_ctr));
        Point2D p = v->get_cp();
        out.write((const char*) &p, sizeof(p));
        DataPoint dp = v->get_data_point();
        data_len = dp->length;
        out.write((const char*) &data_len, sizeof(data_len));
        out.write((const char*) dp->fields, data_len * sizeof(*dp->fields));
        nodes_ctr++;
    }
    //bezier edge nodes
    for(bez_e_it = bezier_mesh->get_edges_begin();
        bez_e_it != bezier_mesh->get_edges_end(); ++bez_e_it){
        const BezierEdge *e = *bez_e_it;
        node_idx[e->get_control_point(1)] = nodes_ctr;
        out.write( (const char*) &nodes_ctr, sizeof(nodes_ctr));
        out.write( (const char*) &(e->get_cp(1).coords), sizeof(Point2D));
        const LinearData& data = e->get_dp(1);
        unsigned data_len = data.length;
        out.write( (const char*) &data_len, sizeof(data_len));
        out.write( (const char*) data.fields, data_len * sizeof(*data.fields));
        nodes_ctr++;
    }
    //cout<<"Wrote "<<nodes_ctr-1<<" nodes to binary file"<<endl;

    /* Write Boundary Vertexs */
    unsigned total_num_bdry_verts = bdry_mesh->get_num_vertices();
    out.write((char*) &total_num_bdry_verts, sizeof(total_num_bdry_verts));

    for(bdry_v_it = bdry_mesh->get_vertices_begin();
        bdry_v_it != bdry_mesh->get_vertices_end(); ++bdry_v_it) {
        BoundaryVertex *bdry_v = *bdry_v_it;
        bdryvert_idx[ bdry_v ] = bdryverts_ctr;

        //write idx
        out.write( (const char*) &bdryverts_ctr, sizeof(bdryverts_ctr));

        //write node idx
        node_num = node_idx[ bdry_v->get_control_point() ];
        out.write( (const char*) &node_num, sizeof(node_num));

        //write fixed
        fixed = bdry_v->get_fixed();
        out.write( (const char*) &fixed, sizeof(fixed));

        bdryverts_ctr++;
    }
    //cout<<"Wrote "<<bdryverts_ctr-1<<" BoundaryVertexs to binary file"<<endl;

    /* Write Boundary Edges */
    unsigned total_num_bdry_edges = bdry_mesh->get_num_edges();
    out.write((char*) &total_num_bdry_edges, sizeof(total_num_bdry_edges));

    for(bdry_e_it = bdry_mesh->get_edges_begin();
        bdry_e_it != bdry_mesh->get_edges_end(); ++bdry_e_it){
        BoundaryEdge *bdry_e = *bdry_e_it;
        QBSpline *spline;
        unsigned num_internal_deboor, num_knots;
        double knot_val;
        double restlength;

        bdryedge_idx[ bdry_e ] = bdryedges_ctr;

        //write bdry edge index
        out.write( (const char*) &bdryedges_ctr, sizeof(unsigned));

        //write first endpoint
        v_num = bdryvert_idx[ bdry_e->get_vertex(0) ];
        out.write( (const char*) &v_num, sizeof(unsigned));

        //write second endpoint
        v_num = bdryvert_idx[ bdry_e->get_vertex(1) ];
        out.write( (const char*) &v_num, sizeof(unsigned));

        //write fixed and color
        spline = bdry_e->get_spline();
        fixed = bdry_e->get_fixed();
        color = bdry_e->get_color();
        restlength = spline->restlength;
        out.write( (const char*) &fixed, sizeof(fixed));
        out.write( (const char*) &color, sizeof(color));
        out.write( (const char*) &restlength, sizeof(restlength));


        //write the internal deboor point node numbers
        //since we already have the endpoints recorded, there is no reason to
        //write them again.  This assumption is violated during the construction
        // of the mesh using the conformal mesher, but we don't care enough to
        // deal with that case: the boundary edges will look like they all have
        // a control point at the origin.
        num_internal_deboor = spline->get_num_deboor()-2;
        out.write( (const char*) &num_internal_deboor, sizeof(unsigned));
        for(i=0; i < num_internal_deboor; i++){
            node_num = node_idx[ spline->get_deboor(i+1) ];
            out.write( (const char*) &node_num, sizeof(unsigned));
        }

        //write num_knots and knot values
        num_knots = spline->k.size();
        out.write( (const char*) &num_knots, sizeof(unsigned));
        for(i=0; i < num_knots; i++){
            knot_val = spline->k[i];
            out.write( (const char*) &knot_val, sizeof(double));
        }

        bdryedges_ctr++;
    }
    //cout<<"Wrote "<<bdryedges_ctr-1<<" BoundaryEdges to binary file"<<endl;


    /* Write Boundary Faces */
    unsigned total_num_bdry_faces = bdry_mesh->get_num_faces();
    out.write((char*) &total_num_bdry_faces, sizeof(total_num_bdry_faces));

    /* Index zero is the null face. */
    bdryface_idx[ NULL ] = 0;
    bdryfaces_ctr = 1;

    for(bdry_f_it = bdry_mesh->get_faces_begin();
        bdry_f_it != bdry_mesh->get_faces_end(); ++bdry_f_it){
        BoundaryFace *bdry_f = *bdry_f_it;
        bdryface_idx[ bdry_f ] = bdryfaces_ctr;

        //write index
        out.write( (const char*) &bdryfaces_ctr, sizeof(unsigned));

        //write min_angle
        double min_angle = bdry_f->get_min_angle();
        out.write( (const char*) &min_angle, sizeof(double));

        //write color
        unsigned color = bdry_f->get_color();
        out.write( (const char*) &color, sizeof(color));

        //write edge idxs
        num_edges = bdry_f->num_edges();
        out.write( (const char*) &num_edges, sizeof(unsigned));
        for(BoundaryFace::edge_iterator it = bdry_f->begin_edges();
            it != bdry_f->end_edges(); ++it) {
          unsigned e = bdryedge_idx[ *it ];
          out.write( (const char*) &e, sizeof(e));
        }

        bdryfaces_ctr++;
    }
    //cout<<"Wrote "<<bdryfaces_ctr-1<<" BoundaryFaces to binary file"<<endl;


    /* Write Bezier Vertexs */
    unsigned total_num_bez_verts = bezier_mesh->get_num_vertices();
    out.write((char*) &total_num_bez_verts, sizeof(total_num_bez_verts));

    for(bez_v_it = bezier_mesh->get_vertices_begin();
        bez_v_it != bezier_mesh->get_vertices_end(); ++bez_v_it){
        BezierVertex *bez_v = *bez_v_it;
        unsigned bdry_type;
        unsigned bdry_idx;
        unsigned bdry_edge_knot_idx;

        bezvert_idx[ bez_v ] = bezverts_ctr;
        //write bezier vertex idx
        out.write( (const char*) &bezverts_ctr, sizeof(unsigned));

        //write node idx
        node_num = node_idx[ bez_v->get_control_point() ];
        out.write( (const char*) &node_num, sizeof(unsigned));

        //write bdry type (0=none, 1=bdry vertex, 2=bdry edge)
        //write bdry idx (depends on bdry_type)
        //write bdry knot idx
        // Who knows what possessed Mark to use that setting of 0/1/2;
        // now I translate from the containment dimension, to maintain
        // binary compatibility.
        switch(bez_v->containment_dimension()) {
          case 0:
            bdry_type = 1;
            bdry_idx= bdryvert_idx[ bez_v->get_bdry_vertex() ];
            bdry_edge_knot_idx=0;
            break;
          case 1:
            bdry_type = 2;
            bdry_idx = bdryedge_idx[ bez_v->get_bdry_edge() ];
            bdry_edge_knot_idx =
              bez_v->get_bdry_edge()->get_spline()->get_bezier_vertex_idx(bez_v);
            assert(bdry_edge_knot_idx > 0); //this index cannot be endpoints
            break;
          case 2:
            bdry_type = 0;
            bdry_idx=0;
            bdry_edge_knot_idx=0;
            break;
          default: assert(0);
        }
        out.write( (const char*) &bdry_type, sizeof(bdry_type));
        out.write( (const char*) &bdry_idx, sizeof(unsigned));
        out.write( (const char*) &bdry_edge_knot_idx, sizeof(unsigned));

        bezverts_ctr++;
    }
    //cout<<"Wrote "<<bezverts_ctr-1<<" BezierVertexs to binary file"<<endl;


    /* Write Bezier Edges */
    unsigned total_num_bez_edges = bezier_mesh->get_num_edges();
    out.write((char*) &total_num_bez_edges, sizeof(total_num_bez_edges));

    for(bez_e_it = bezier_mesh->get_edges_begin();
        bez_e_it != bezier_mesh->get_edges_end(); ++bez_e_it){
        BezierEdge *bez_e = *bez_e_it;
        unsigned bdry_e_num, bdry_u0_idx, bdry_u1_idx;
        BezierVertex *v0, *v1;
        bool orientation;
        unsigned last_idx;
        bezedge_idx[ bez_e ] = bezedges_ctr;

        //write bezier edge idx
        out.write( (const char*) &bezedges_ctr, sizeof(unsigned));

        //write endpoint idxs
        v_num = bezvert_idx[ bez_e->get_vertex(0) ];
        out.write( (const char*) &v_num, sizeof(unsigned));
        v_num = bezvert_idx[ bez_e->get_vertex(1) ];
        out.write( (const char*) &v_num, sizeof(unsigned));

        //write node idx
        node_num = node_idx[ bez_e->get_control_point(1) ];
        out.write( (const char*) &node_num, sizeof(unsigned));

        //write bdry edge idx, (0 if not on bdry)
        if(bez_e->is_boundary()) {
            BoundaryEdge *bdry_e = bez_e->get_bdry_edge();
            bdry_e_num = bdryedge_idx[ bdry_e ];
            v0 = bez_e->get_vertex(0);
            v1 = bez_e->get_vertex(1);

            /* This will get the knot index if the vertex is a D1*/
            bdry_u0_idx = bdry_e->get_spline()->get_bezier_vertex_idx( v0 );
            bdry_u1_idx = bdry_e->get_spline()->get_bezier_vertex_idx( v1 );

            /* Otherwise the vertex is a D0 vertex (a endpoint of the spline)
             * so find it's knot index with respect to the endpoints of the Boundary Edge
             * again we avoid comparison of floating-points that would otherwise be necessary
             */

            /* determine how spline is parameterized with respecet to the endpoints of the bdry edge
             * if the first vertex of the bdry edge cooresponds to the begining of the parameterization
             * then orientation is true
             */
            QBSpline *spline = bdry_e->get_spline();
            orientation = (spline->b[0] == bdry_e->get_vertex(0)->get_control_point());
            last_idx = spline->k.size()-1;

            //just a check.  if both of bezier edge's vertexs are endpoints of spline, then spline must be short
            if(bdry_u0_idx ==0 && bdry_u1_idx==0){
                assert(spline->k.size()==2);
                assert(!spline->closed);
            }

            //if bezier edge's first vertex is endpoint, then find it's knot idx
            if(bdry_u0_idx == 0){
                if(spline->closed){
                    if(bdry_u1_idx-1 == 0) bdry_u0_idx=0;
                    else bdry_u0_idx=last_idx;
                } else {
                    if( v0 == bdry_e->get_vertex(0)->get_bezier())
                        bdry_u0_idx = (orientation) ? 0 : last_idx;
                    else if( v0 == bdry_e->get_vertex(1)->get_bezier())
                        bdry_u0_idx = (orientation) ? last_idx : 0;
                }
            }

            //if bezier edge's second vertex is endpoint, then find it's knot idx
            if(bdry_u1_idx == 0){
                if(spline->closed){
                    if(bdry_u0_idx-1 == 0) bdry_u1_idx=0;
                    else bdry_u1_idx=last_idx;
                } else{
                    if( v1 == bdry_e->get_vertex(0)->get_bezier())
                        bdry_u1_idx = (orientation) ? 0 : last_idx;
                    else if( v1 == bdry_e->get_vertex(1)->get_bezier())
                        bdry_u1_idx = (orientation) ? last_idx : 0;
                }
            }
            //cout<<"orientation: "<<orientation<<" u0_idx:"<<bdry_u0_idx<<" u1_idx:"<<bdry_u1_idx<<" last_idx:"<<last_idx<<endl;
            assert( abs((int)bdry_u0_idx - (int)bdry_u1_idx) == 1);

        } else {
            bdry_e_num = 0;
            bdry_u0_idx = 0;
            bdry_u1_idx = 0;
        }

        out.write( (const char*) &bdry_e_num,  sizeof(unsigned));
        out.write( (const char*) &bdry_u0_idx, sizeof(unsigned));
        out.write( (const char*) &bdry_u1_idx, sizeof(unsigned));

        bezedges_ctr++;
    }
    //cout<<"[MeshIO.C] Wrote "<<bezedges_ctr-1<<" BezierEdges to binary file"<<endl;

    /* Write Bezier Triangles */
    unsigned total_num_bez_tris = bezier_mesh->get_num_faces();
    out.write((char*) &total_num_bez_tris, sizeof(total_num_bez_tris));


    //cout << "[MeshIO.C] Writing bezier triangles"<<endl;
    for(bez_t_it = bezier_mesh->get_faces_begin();
        bez_t_it != bezier_mesh->get_faces_end(); ++bez_t_it){
        BezierTriangle *bez_t = *bez_t_it;

        //write edge idxs
        e_num = bezedge_idx[ bez_t->get_edge(0) ];
        out.write( (const char*) &e_num, sizeof(unsigned));
        e_num = bezedge_idx[ bez_t->get_edge(1) ];
        out.write( (const char*) &e_num, sizeof(unsigned));
        e_num = bezedge_idx[ bez_t->get_edge(2) ];
        out.write( (const char*) &e_num, sizeof(unsigned));

    //write inverse_handed_flag, tp517:4/16/07
    bool ihf = bezier_mesh->is_inverse_handed(bez_t);
    out.write( (const char*) &ihf, sizeof(bool));

        //write bdry face
        f_num = bdryface_idx[ bez_t->get_bdry_face_or_null() ];
        out.write( (const char*) &f_num, sizeof(unsigned));

        beztris_ctr++;
    }
    //cout<<"Wrote "<<beztris_ctr-1<<" BezierTriangles to binary file"<<endl;

    out.close();

    return true;
}

bool MeshBinaryOutput::write(const char *filename, BezierMesh *bez,
                             BoundaryMesh *bdry)
{
  MeshBinaryOutput mo(bez,bdry);
  return mo.write(filename);
}

bool MeshBinaryOutput::write(const char *filename, Simulation *sim)
{
  MeshBinaryOutput mo(sim);
  return mo.write(filename);
}

---