MOAB: Mesh Oriented datABase  (version 5.2.1)
MeshImpl.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government 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     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov,
00025     kraftche@cae.wisc.edu
00026 
00027   ***************************************************************** */
00028 //
00029 //   SUMMARY:
00030 //     USAGE:
00031 //
00032 // ORIG-DATE: 16-May-02 at 10:26:21
00033 //  LAST-MOD: 15-Nov-04 by kraftche@cae.wisc.edu
00034 //
00035 /*! \file MeshImpl.cpp
00036 
00037 \brief This files contains a mesh database implementation that can be used
00038 to run mesquite by default.
00039 
00040     \author Thomas Leurent
00041     \author Darryl Melander
00042     \author Jason Kraftcheck
00043     \date 2004-11-15
00044  */
00045 
00046 #include "MeshImpl.hpp"
00047 #include "FileTokenizer.hpp"
00048 #include "Vector3D.hpp"
00049 #include "MsqVertex.hpp"
00050 #include "MeshImplData.hpp"
00051 #include "MeshImplTags.hpp"
00052 #include "MsqDebug.hpp"
00053 #include "MsqError.hpp"
00054 #include "VtkTypeInfo.hpp"
00055 
00056 #include <string>
00057 #include <vector>
00058 #include <algorithm>
00059 #include <functional>
00060 #include <map>
00061 using std::string;
00062 using std::vector;
00063 
00064 #include <fstream>
00065 #include <iomanip>
00066 using std::endl;
00067 using std::ifstream;
00068 using std::ofstream;
00069 
00070 #ifdef MSQ_USING_EXODUS
00071 #include "exodusII.h"
00072 #endif
00073 
00074 #include "MsqDebug.hpp"
00075 namespace MBMesquite
00076 {
00077 
00078 const char* MESQUITE_FIELD_TAG = "MesquiteTags";
00079 
00080 MeshImpl::MeshImpl() : numCoords( 0 ), myMesh( new MeshImplData ), myTags( new MeshImplTags ) {}
00081 
00082 MeshImpl::~MeshImpl()
00083 {
00084     delete myMesh;
00085     delete myTags;
00086 }
00087 
00088 MeshImpl::MeshImpl( int num_nodes, int num_elem, EntityTopology entity_topology, const bool* fixed,
00089                     const double* coords, const int* connectivity )
00090     : numCoords( 3 ), myMesh( new MeshImplData ), myTags( new MeshImplTags )
00091 {
00092     MsqError err;
00093     myMesh->allocate_vertices( num_nodes, err );MSQ_ERRRTN( err );
00094     myMesh->allocate_elements( num_elem, err );MSQ_ERRRTN( err );
00095 
00096     // Fill in the data
00097     if( fixed )
00098     {
00099         for( int i = 0; i < num_nodes; ++i )
00100         {
00101             myMesh->reset_vertex( i, Vector3D( coords + 3 * i ), fixed[i], err );
00102         }
00103     }
00104     else
00105     {
00106         for( int i = 0; i < num_nodes; ++i )
00107         {
00108             myMesh->reset_vertex( i, Vector3D( coords + 3 * i ), false, err );
00109         }
00110     }
00111 
00112     const int verts_per_elem = TopologyInfo::corners( entity_topology );
00113 
00114     std::vector< long > connect( verts_per_elem );
00115     const int* conn_iter = connectivity;
00116     for( int i = 0; i < num_elem; i++, conn_iter += verts_per_elem )
00117     {
00118         std::copy( conn_iter, conn_iter + verts_per_elem, connect.begin() );
00119         myMesh->reset_element( i, connect, entity_topology, err );
00120     }
00121 }
00122 
00123 MeshImpl::MeshImpl( int num_nodes, int num_elem, const EntityTopology* element_topologies, const bool* fixed,
00124                     const double* coords, const int* connectivity )
00125     : numCoords( 3 ), myMesh( new MeshImplData ), myTags( new MeshImplTags )
00126 {
00127     MsqError err;
00128     myMesh->allocate_vertices( num_nodes, err );MSQ_ERRRTN( err );
00129     myMesh->allocate_elements( num_elem, err );MSQ_ERRRTN( err );
00130 
00131     // Fill in the data
00132     if( fixed )
00133     {
00134         for( int i = 0; i < num_nodes; ++i )
00135         {
00136             myMesh->reset_vertex( i, Vector3D( coords + 3 * i ), fixed[i], err );
00137         }
00138     }
00139     else
00140     {
00141         for( int i = 0; i < num_nodes; ++i )
00142         {
00143             myMesh->reset_vertex( i, Vector3D( coords + 3 * i ), false, err );
00144         }
00145     }
00146 
00147     int num_indices = 0;
00148     std::vector< long > connect;
00149 
00150     // Count the number of indices
00151     const int* conn_iter = connectivity;
00152     for( int i = 0; i < num_elem; ++i )
00153     {
00154         num_indices = TopologyInfo::corners( element_topologies[i] );
00155         connect.resize( num_indices );
00156         std::copy( conn_iter, conn_iter + num_indices, connect.begin() );
00157         myMesh->reset_element( i, connect, element_topologies[i], err );
00158         conn_iter += num_indices;
00159     }
00160 }
00161 
00162 void MeshImpl::clear()
00163 {
00164     myMesh->clear();
00165     myTags->clear();
00166 }
00167 
00168 void MeshImpl::set_all_fixed_flags( bool value, MsqError& err )
00169 {
00170     for( size_t i = 0; i < myMesh->max_vertex_index(); ++i )
00171     {
00172         if( myMesh->is_vertex_valid( i ) )
00173         {
00174             myMesh->fix_vertex( i, value, err );MSQ_ERRRTN( err );
00175         }
00176     }
00177 }
00178 
00179 void MeshImpl::set_all_slaved_flags( bool, MsqError& err )
00180 {
00181     for( size_t e = 0; e < myMesh->max_element_index(); ++e )
00182     {
00183         if( !myMesh->is_element_valid( e ) ) continue;
00184 
00185         // Get element connectivity
00186         const std::vector< size_t >& verts = myMesh->element_connectivity( e, err );MSQ_ERRRTN( err );
00187 
00188         // Get element properties
00189         EntityTopology type = myMesh->element_topology( e, err );MSQ_ERRRTN( err );
00190         unsigned ncorner = TopologyInfo::corners( type );
00191 
00192         for( unsigned i = 0; i < ncorner; ++i )
00193         {
00194             myMesh->slave_vertex( verts[i], false, err );MSQ_ERRRTN( err );
00195         }
00196         for( unsigned i = ncorner; i < verts.size(); ++i )
00197         {
00198             myMesh->slave_vertex( verts[i], false, err );MSQ_ERRRTN( err );
00199         }
00200     }
00201 }
00202 
00203 /**\brief Helper function for MeshImpl::mark_skin_fixed */
00204 static bool is_side_boundary( MeshImplData* myMesh, size_t elem, unsigned side_dim, unsigned side_num, MsqError& err )
00205 {
00206     // Get the vertices of the side as indices into the above 'verts' list.
00207     const EntityTopology type = myMesh->element_topology( elem, err );
00208     MSQ_ERRZERO( err );
00209     unsigned n;  // number of vertices
00210     const unsigned* conn = TopologyInfo::side_vertices( type, side_dim, side_num, n, err );
00211     MSQ_ERRZERO( err );
00212 
00213     // start with the assumption that the side is on the bounary
00214     bool boundary = true;
00215 
00216     // get vertices in element connectivity
00217     const std::vector< size_t >& verts = myMesh->element_connectivity( elem, err );
00218     MSQ_ERRZERO( err );
00219 
00220     // Choose one vertex in face, and get adjacent elements to that vertex
00221     const std::vector< size_t >& elems = myMesh->vertex_adjacencies( verts[conn[0]], err );
00222     MSQ_ERRZERO( err );
00223 
00224     // For each adjacent element
00225     for( unsigned i = 0; i < elems.size(); ++i )
00226     {
00227         // we want *other* adjacent elements
00228         if( elems[i] == elem ) continue;
00229 
00230         // skip elements of smaller dimension
00231         EntityTopology type2 = myMesh->element_topology( elems[i], err );
00232         if( TopologyInfo::dimension( type2 ) <= side_dim ) continue;
00233 
00234         // get number of 'sides' of the appropriate dimension.
00235         const std::vector< size_t >& verts2 = myMesh->element_connectivity( elems[i], err );
00236         MSQ_ERRZERO( err );
00237         int sides2 = TopologyInfo::adjacent( type2, side_dim );
00238         for( int j = 0; j < sides2; ++j )
00239         {
00240             if( TopologyInfo::compare_sides( (const size_t*)arrptr( verts ), type, side_num,
00241                                              (const size_t*)arrptr( verts2 ), type2, j, side_dim, err ) )
00242                 boundary = false;
00243             MSQ_ERRZERO( err );
00244         }
00245     }
00246 
00247     return boundary;
00248 }
00249 
00250 void MeshImpl::set_skin_flags( bool corner_fixed_flag, bool midnode_fixed_flag, bool midnode_slaved_flag,
00251                                MsqError& err )
00252 {
00253     // For each element, for each side of that element, check for
00254     // an adjacent element.
00255     for( size_t i = 0; i < myMesh->max_element_index(); ++i )
00256     {
00257         if( !myMesh->is_element_valid( i ) ) continue;
00258 
00259         // Get element connectivity
00260         const std::vector< size_t >& verts = myMesh->element_connectivity( i, err );
00261 
00262         // Get element properties
00263         EntityTopology type = myMesh->element_topology( i, err );
00264         unsigned dim        = TopologyInfo::dimension( type );
00265         int sides           = TopologyInfo::adjacent( type, dim - 1 );
00266         bool midedge, midface, midvol;
00267         TopologyInfo::higher_order( type, verts.size(), midedge, midface, midvol, err );MSQ_ERRRTN( err );
00268         const bool midside    = ( dim == 2 && midedge ) || ( dim == 3 && midface );
00269         const bool midsubside = dim == 3 && midedge;
00270 
00271         // For each side of the element (each edge for surface elems,
00272         // each face for volume elements)..
00273         for( int j = 0; j < sides; ++j )
00274         {
00275             // Get the vertices of the side as indices into the above 'verts' list.
00276             unsigned n;  // number of vertices
00277             const unsigned* conn = TopologyInfo::side_vertices( type, dim - 1, j, n, err );MSQ_ERRRTN( err );
00278 
00279             // if side is on boundary, mark side vertices appropriately
00280             bool boundary = is_side_boundary( myMesh, i, dim - 1, j, err );MSQ_ERRRTN( err );
00281             if( boundary )
00282             {
00283                 // mark corner vertices as fixed
00284                 for( unsigned k = 0; k < n; ++k )
00285                 {
00286                     myMesh->fix_vertex( verts[conn[k]], corner_fixed_flag, err );MSQ_ERRRTN( err );
00287                 }
00288 
00289                 // mark higher-order node in center of side as fixed
00290                 if( midside )
00291                 {
00292                     unsigned idx = TopologyInfo::higher_order_from_side( type, verts.size(), dim - 1, j, err );MSQ_ERRRTN( err );
00293                     myMesh->fix_vertex( verts[idx], midnode_fixed_flag, err );MSQ_ERRRTN( err );
00294                     myMesh->slave_vertex( verts[idx], midnode_slaved_flag, err );MSQ_ERRRTN( err );
00295                 }
00296 
00297                 // if side is a face, mark nodes on edges of face as fixed
00298                 if( midsubside )
00299                 {
00300                     for( unsigned k = 0; k < n; ++k )
00301                     {
00302                         unsigned edge[2] = { conn[k], conn[( k + 1 ) % n] };
00303                         bool r;
00304                         unsigned edge_num = TopologyInfo::find_edge( type, edge, r, err );MSQ_ERRRTN( err );
00305 
00306                         unsigned idx = TopologyInfo::higher_order_from_side( type, verts.size(), 1, edge_num, err );MSQ_ERRRTN( err );
00307                         myMesh->fix_vertex( verts[idx], midnode_fixed_flag, err );MSQ_ERRRTN( err );
00308                         myMesh->slave_vertex( verts[idx], midnode_slaved_flag, err );MSQ_ERRRTN( err );
00309                     }
00310                 }
00311             }
00312         }  // for (j in sides)
00313     }      // for (i in elems)
00314 }
00315 
00316 void MeshImpl::mark_skin_fixed( MsqError& err, bool clear_existing )
00317 {
00318     if( clear_existing )
00319     {
00320         set_all_fixed_flags( false, err );MSQ_ERRRTN( err );
00321     }
00322 
00323     set_skin_flags( true, true, false, err );MSQ_ERRRTN( err );
00324 }
00325 
00326 static void get_field_names( const TagDescription& tag, std::string& field_out, std::string& member_out, MsqError& err )
00327 {
00328     std::string::size_type idx;
00329 
00330     if( tag.vtkType != TagDescription::FIELD )
00331     {
00332         field_out  = tag.name;
00333         member_out = "";
00334     }
00335     else if( !tag.member.empty() )
00336     {
00337         field_out  = tag.name;
00338         member_out = tag.member;
00339     }
00340     // If field is a struct, then the Mesquite tag name is a contcatenation
00341     // of the field name and the struct member name.  Extract them.
00342     else
00343     {
00344         idx = tag.name.find( " " );
00345         if( idx == std::string::npos )
00346         {
00347             field_out  = tag.name;
00348             member_out = MESQUITE_FIELD_TAG;
00349         }
00350         else
00351         {
00352             field_out  = tag.name.substr( 0, idx );
00353             member_out = tag.name.substr( idx + 1 );
00354         }
00355     }
00356 
00357     idx = field_out.find( " " );
00358     if( idx != std::string::npos ) MSQ_SETERR( err )
00359     ( MsqError::FILE_FORMAT, "Cannot write tag name \"%s\" containing spaces to VTK file.", field_out.c_str() );
00360 
00361     idx = member_out.find( " " );
00362     if( idx != std::string::npos ) MSQ_SETERR( err )
00363     ( MsqError::FILE_FORMAT, "Cannot write field member name \"%s\" containing spaces to VTK file.",
00364       member_out.c_str() );
00365 }
00366 
00367 void MeshImpl::write_vtk( const char* out_filename, MsqError& err )
00368 {
00369     ofstream file( out_filename );
00370     if( !file )
00371     {
00372         MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00373         return;
00374     }
00375 
00376     file.precision( 10 );
00377 
00378     // Write a header
00379     file << "# vtk DataFile Version 3.0\n";
00380     file << "Mesquite Mesh\n";
00381     file << "ASCII\n";
00382     file << "DATASET UNSTRUCTURED_GRID\n";
00383 
00384     // Write vertex coordinates
00385     file << "POINTS " << myMesh->num_vertices() << " double\n";
00386 
00387     std::vector< size_t > vertex_indices( myMesh->max_vertex_index() );
00388     size_t i, count = 0;
00389     for( i = 0; i < myMesh->max_vertex_index(); ++i )
00390     {
00391         if( myMesh->is_vertex_valid( i ) )
00392         {
00393             Vector3D coords = myMesh->get_vertex_coords( i, err );MSQ_ERRRTN( err );
00394             file << coords[0] << ' ' << coords[1] << ' ' << coords[2] << '\n';
00395             vertex_indices[i] = count++;
00396         }
00397         else
00398         {
00399             vertex_indices[i] = myMesh->max_vertex_index();
00400         }
00401     }
00402 
00403     // Write out the connectivity table
00404     size_t elem_idx;
00405     size_t connectivity_size = myMesh->num_elements() + myMesh->num_vertex_uses();
00406     file << "CELLS " << myMesh->num_elements() << ' ' << connectivity_size << '\n';
00407     for( elem_idx = 0; elem_idx < myMesh->max_element_index(); ++elem_idx )
00408     {
00409         if( !myMesh->is_element_valid( elem_idx ) ) continue;
00410 
00411         std::vector< size_t > conn = myMesh->element_connectivity( elem_idx, err );MSQ_ERRRTN( err );
00412         EntityTopology topo = myMesh->element_topology( elem_idx, err );MSQ_ERRRTN( err );
00413 
00414         // If necessary, convert from Exodus to VTK node-ordering.
00415         const VtkTypeInfo* info = VtkTypeInfo::find_type( topo, conn.size(), err );MSQ_ERRRTN( err );
00416         if( info->msqType != POLYGON ) info->mesquiteToVtkOrder( conn );
00417 
00418         file << conn.size();
00419         for( i = 0; i < conn.size(); ++i )
00420             file << ' ' << vertex_indices[(size_t)conn[i]];
00421         file << '\n';
00422     }
00423 
00424     // Write out the element types
00425     file << "CELL_TYPES " << myMesh->num_elements() << '\n';
00426     for( elem_idx = 0; elem_idx < myMesh->max_element_index(); ++elem_idx )
00427     {
00428         if( !myMesh->is_element_valid( elem_idx ) ) continue;
00429 
00430         EntityTopology topo = myMesh->element_topology( elem_idx, err );MSQ_ERRRTN( err );
00431         count = myMesh->element_connectivity( elem_idx, err ).size();MSQ_ERRRTN( err );
00432         const VtkTypeInfo* info = VtkTypeInfo::find_type( topo, count, err );MSQ_ERRRTN( err );
00433         file << info->vtkType << '\n';
00434     }
00435 
00436     // Write out which points are fixed.
00437     file << "POINT_DATA " << myMesh->num_vertices() << "\nSCALARS fixed int\nLOOKUP_TABLE default\n";
00438     for( i = 0; i < myMesh->max_vertex_index(); ++i )
00439         if( myMesh->is_vertex_valid( i ) ) file << ( myMesh->vertex_is_fixed( i, err ) ? "1" : "0" ) << "\n";
00440 
00441     if( myMesh->have_slaved_flags() )
00442     {
00443         file << "SCALARS slaved int\nLOOKUP_TABLE default\n";
00444         for( i = 0; i < myMesh->max_vertex_index(); ++i )
00445             if( myMesh->is_vertex_valid( i ) ) file << ( myMesh->vertex_is_slaved( i, err ) ? "1" : "0" ) << "\n";
00446     }
00447 
00448     // Make pass over the list of tags to:
00449     // - Check if there are any tags on elements.
00450     // - Get the list of field names by which to group tag data.
00451     // - Assign VTK types for tags that have an unknown VTK type.
00452     MeshImplTags::TagIterator tagiter;
00453     std::multimap< std::string, size_t > fields;
00454     bool have_elem_data = false;
00455     for( tagiter = myTags->tag_begin(); tagiter != myTags->tag_end(); ++tagiter )
00456     {
00457         bool haveelem = myTags->tag_has_element_data( *tagiter, err );MSQ_ERRRTN( err );
00458         if( haveelem ) have_elem_data = true;
00459 
00460         const TagDescription& desc = myTags->properties( *tagiter, err );
00461         std::string field, member;
00462         get_field_names( desc, field, member, err );MSQ_ERRRTN( err );
00463         fields.insert( std::multimap< std::string, size_t >::value_type( field, *tagiter ) );
00464     }
00465 
00466     // Write vertex tag data to vtk attributes, group by field.
00467     // We already wrote the header line for the POINT_DATA block
00468     // with the fixed flag, so just write the rest of the tags now.
00469     std::multimap< std::string, size_t >::iterator f, e;
00470     f = fields.begin();
00471     while( f != fields.end() )
00472     {
00473         if( !myTags->tag_has_vertex_data( f->second, err ) )
00474         {
00475             ++f;
00476             continue;
00477         }
00478 
00479         int pcount = 0;
00480         for( e = f; e != fields.end() && e->first == f->first; ++e )
00481             if( myTags->tag_has_vertex_data( e->second, err ) ) ++pcount;
00482         if( !pcount ) continue;
00483 
00484         if( myTags->properties( f->second, err ).vtkType == TagDescription::FIELD )
00485             file << "FIELD " << f->first << " " << pcount << std::endl;
00486         else if( pcount > 1 )
00487         {
00488             MSQ_SETERR( err )
00489             ( MsqError::INTERNAL_ERROR, "Tag name \"%s\" conflicts with VTK field name in tag \"%s\"\n",
00490               f->first.c_str(), myTags->properties( ( ++f )->second, err ).name.c_str() );
00491             return;
00492         }
00493 
00494         for( ; f != e; ++f )
00495         {
00496             if( !myTags->tag_has_vertex_data( f->second, err ) ) continue;
00497 
00498             const TagDescription& desc = myTags->properties( f->second, err );MSQ_ERRRTN( err );
00499             std::vector< char > tagdata( myMesh->num_vertices() * desc.size );
00500             std::vector< char >::iterator iter = tagdata.begin();
00501             for( i = 0; i < myMesh->max_vertex_index(); ++i )
00502             {
00503                 if( myMesh->is_vertex_valid( i ) )
00504                 {
00505                     myTags->get_vertex_data( f->second, 1, &i, &*iter, err );
00506                     if( err.error_code() == MsqError::TAG_NOT_FOUND )
00507                     {
00508                         memset( &*iter, 0, desc.size );
00509                         err.clear();
00510                     }
00511                     else if( MSQ_CHKERR( err ) )
00512                         return;
00513                     iter += desc.size;
00514                 }
00515             }
00516             vtk_write_attrib_data( file, desc, arrptr( tagdata ), myMesh->num_vertices(), err );MSQ_ERRRTN( err );
00517         }
00518     }
00519 
00520     // If no element tags, then done
00521     if( !have_elem_data )
00522     {
00523         file.close();
00524         return;
00525     }
00526 
00527     // Begin element data section
00528     file << "\nCELL_DATA " << myMesh->num_elements() << "\n";
00529     // Write element tag data to vtk attributes, group by field.
00530     f = fields.begin();
00531     while( f != fields.end() )
00532     {
00533         if( !myTags->tag_has_element_data( f->second, err ) )
00534         {
00535             ++f;
00536             continue;
00537         }
00538 
00539         int pcount = 0;
00540         for( e = f; e != fields.end() && e->first == f->first; ++e )
00541             if( myTags->tag_has_element_data( e->second, err ) ) ++pcount;
00542         if( !pcount ) continue;
00543 
00544         if( myTags->properties( f->second, err ).vtkType == TagDescription::FIELD )
00545             file << "FIELD " << f->first << " " << pcount << std::endl;
00546         else if( pcount > 1 )
00547         {
00548             MSQ_SETERR( err )
00549             ( MsqError::INTERNAL_ERROR, "Tag name \"%s\" conflicts with VTK field name in tag \"%s\"\n",
00550               f->first.c_str(), myTags->properties( ( ++f )->second, err ).name.c_str() );
00551             return;
00552         }
00553 
00554         for( ; f != e; ++f )
00555         {
00556             if( !myTags->tag_has_element_data( f->second, err ) ) continue;
00557 
00558             const TagDescription& desc = myTags->properties( f->second, err );MSQ_ERRRTN( err );
00559             std::vector< char > tagdata( myMesh->num_elements() * desc.size );
00560             std::vector< char >::iterator iter = tagdata.begin();
00561             for( i = 0; i < myMesh->max_element_index(); ++i )
00562             {
00563                 if( myMesh->is_element_valid( i ) )
00564                 {
00565                     myTags->get_element_data( f->second, 1, &i, &*iter, err );
00566                     if( err.error_code() == MsqError::TAG_NOT_FOUND )
00567                     {
00568                         memset( &*iter, 0, desc.size );
00569                         err.clear();
00570                     }
00571                     else if( MSQ_CHKERR( err ) )
00572                         return;
00573                     iter += desc.size;
00574                 }
00575             }
00576             vtk_write_attrib_data( file, desc, arrptr( tagdata ), myMesh->num_elements(), err );MSQ_ERRRTN( err );
00577         }
00578     }
00579 
00580     // Close the file
00581     file.close();
00582 }
00583 
00584 void MeshImpl::read_exodus( const char* in_filename, MsqError& err )
00585 {
00586 #ifndef MSQ_USING_EXODUS
00587     MSQ_SETERR( err )( MsqError::NOT_IMPLEMENTED );
00588     MSQ_DBGOUT( 1 ) << "Cannot read ExodusII file: " << in_filename << "\n";
00589     return;
00590 #else
00591 
00592     clear();
00593 
00594     int app_float_size  = sizeof( double );
00595     int file_float_size = 0;
00596     float exo_version   = 0;
00597     int exo_err         = 0;
00598 
00599     // Open the file
00600     int file_id = ex_open( in_filename, EX_READ, &app_float_size, &file_float_size, &exo_version );
00601 
00602     // Make sure we opened the file correctly
00603     if( file_id < 0 )
00604     {
00605         MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00606         return;
00607     }
00608 
00609     // make sure the file is saved as doubles
00610     if( file_float_size != sizeof( double ) )
00611     {
00612         MSQ_SETERR( err )
00613         ( "File saved with float-sized reals.  Can only read files "
00614           "saved with doubles.",
00615           MsqError::NOT_IMPLEMENTED );
00616         return;
00617     }
00618 
00619     char title[MAX_LINE_LENGTH];
00620     int dim, vert_count, elem_count, block_count, ns_count, ss_count;
00621 
00622     // get info about the file
00623     exo_err = ex_get_init( file_id, title, &dim, &vert_count, &elem_count, &block_count, &ns_count, &ss_count );
00624     if( exo_err < 0 )
00625     {
00626         MSQ_SETERR( err )( "Unable to get entity counts from file.", MsqError::PARSE_ERROR );
00627         return;
00628     }
00629 
00630     myMesh->allocate_vertices( vert_count, err );MSQ_ERRRTN( err );
00631     myMesh->allocate_elements( elem_count, err );MSQ_ERRRTN( err );
00632 
00633     // Now fill in the data
00634 
00635     // Get the vertex coordinates
00636     std::vector< double > coords( vert_count * 3 );
00637     double* x_iter = arrptr( coords );
00638     double* y_iter = &coords[vert_count];
00639     double* z_iter = &coords[2 * vert_count];
00640     numCoords      = dim;
00641     if( dim == 2 )
00642     {
00643         exo_err = ex_get_coord( file_id, x_iter, y_iter, 0 );
00644         memset( z_iter, 0, sizeof( double ) * vert_count );
00645     }
00646     else
00647     {
00648         exo_err = ex_get_coord( file_id, x_iter, y_iter, z_iter );
00649     }
00650     // Make sure it worked
00651     if( exo_err < 0 )
00652     {
00653         MSQ_SETERR( err )
00654         ( "Unable to retrieve vertex coordinates from file.", MsqError::PARSE_ERROR );
00655         return;
00656     }
00657 
00658     // Store vertex coordinates in vertex array
00659     int i;
00660     for( i = 0; i < vert_count; ++i )
00661         myMesh->reset_vertex( i, Vector3D( *( x_iter++ ), *( y_iter++ ), *( z_iter )++ ), false, err );
00662     coords.clear();
00663 
00664     // Get block list
00665     std::vector< int > block_ids( block_count );
00666     exo_err = ex_get_elem_blk_ids( file_id, arrptr( block_ids ) );
00667     if( exo_err < 0 )
00668     {
00669         MSQ_SETERR( err )( "Unable to read block IDs from file.", MsqError::PARSE_ERROR );
00670         return;
00671     }
00672 
00673     std::vector< int > conn;
00674     size_t index = 0;
00675     for( i = 0; i < block_count; i++ )
00676     {
00677         // Get info about this block's elements
00678         char elem_type_str[MAX_STR_LENGTH];
00679         int num_block_elems, verts_per_elem, num_atts;
00680         exo_err =
00681             ex_get_elem_block( file_id, block_ids[i], elem_type_str, &num_block_elems, &verts_per_elem, &num_atts );
00682         if( exo_err < 0 )
00683         {
00684             MSQ_SETERR( err )( "Unable to read parameters for block.", MsqError::PARSE_ERROR );
00685             return;
00686         }
00687 
00688         // Figure out which type of element we're working with
00689         EntityTopology elem_type;
00690         for( int j = 0; elem_type_str[j]; j++ )
00691             elem_type_str[j] = toupper( elem_type_str[j] );
00692         if( !strncmp( elem_type_str, "TRI", 3 ) ) { elem_type = TRIANGLE; }
00693         else if( !strncmp( elem_type_str, "QUA", 3 ) || !strncmp( elem_type_str, "SHE", 3 ) )
00694         {
00695             elem_type = QUADRILATERAL;
00696         }
00697         else if( !strncmp( elem_type_str, "HEX", 3 ) )
00698         {
00699             elem_type = HEXAHEDRON;
00700         }
00701         else if( !strncmp( elem_type_str, "TET", 3 ) )
00702         {
00703             elem_type = TETRAHEDRON;
00704         }
00705         else if( !strncmp( elem_type_str, "PYRAMID", 7 ) )
00706         {
00707             elem_type = PYRAMID;
00708         }
00709         else if( !strncmp( elem_type_str, "WEDGE", 5 ) )
00710         {
00711             elem_type = PRISM;
00712         }
00713         else
00714         {
00715             MSQ_SETERR( err )
00716             ( "Unrecognized element type in block", MsqError::UNSUPPORTED_ELEMENT );
00717             continue;
00718         }
00719 
00720         if( conn.size() < (unsigned)( num_block_elems * verts_per_elem ) )
00721             conn.resize( num_block_elems * verts_per_elem );
00722         exo_err = ex_get_elem_conn( file_id, block_ids[i], arrptr( conn ) );
00723         if( exo_err < 0 )
00724         {
00725             MSQ_SETERR( err )
00726             ( "Unable to read element block connectivity.", MsqError::PARSE_ERROR );
00727             return;
00728         }
00729 
00730         std::vector< size_t > vertices( verts_per_elem );
00731         std::vector< int >::iterator conn_iter = conn.begin();
00732         for( const size_t end = index + num_block_elems; index < end; ++index )
00733         {
00734             for( std::vector< size_t >::iterator iter = vertices.begin(); iter != vertices.end(); ++iter, ++conn_iter )
00735                 *iter = *conn_iter - 1;
00736 
00737             myMesh->reset_element( index, vertices, elem_type, err );MSQ_CHKERR( err );
00738         }
00739     }
00740 
00741     // Finally, mark boundary nodes
00742     int num_fixed_nodes = 0;
00743     int num_dist_in_set = 0;
00744     if( ns_count > 0 )
00745     {
00746         exo_err = ex_get_node_set_param( file_id, 111, &num_fixed_nodes, &num_dist_in_set );
00747         if( exo_err < 0 )
00748         {
00749             MSQ_PRINT( 1 )( "\nError opening nodeset 111, no boundary nodes marked." );
00750             num_fixed_nodes = 0;
00751         }
00752     }
00753     std::vector< int > fixed_nodes( num_fixed_nodes );
00754     if( num_fixed_nodes )
00755     {
00756         exo_err = ex_get_node_set( file_id, 111, arrptr( fixed_nodes ) );
00757         if( exo_err < 0 ) { MSQ_SETERR( err )( "Error retrieving fixed nodes.", MsqError::PARSE_ERROR ); }
00758     }
00759 
00760     // See if this vertex is marked as a boundary vertex
00761     for( i = 0; i < num_fixed_nodes; ++i )
00762     {
00763         myMesh->fix_vertex( fixed_nodes[i] - 1, true, err );MSQ_CHKERR( err );
00764     }
00765 
00766     // Finish up
00767     exo_err = ex_close( file_id );
00768     if( exo_err < 0 ) MSQ_SETERR( err )( "Error closing Exodus file.", MsqError::IO_ERROR );
00769 #endif
00770 }
00771 //! Writes an exodus file of the mesh.
00772 void MeshImpl::write_exodus( const char* out_filename, MsqError& err )
00773 {
00774     // just return an error if we don't have access to exodus
00775 #ifndef MSQ_USING_EXODUS
00776     MSQ_SETERR( err )( "Exodus not enabled in this build of Mesquite", MsqError::NOT_IMPLEMENTED );
00777     MSQ_DBGOUT( 1 ) << "Cannot read ExodusII file: " << out_filename << "\n";
00778     return;
00779 #else
00780     size_t i, j, k;
00781     if( !myMesh || !myMesh->num_vertices() )
00782     {
00783         MSQ_SETERR( err )
00784         ( "No vertices in MeshImpl.  Nothing written to file.", MsqError::PARSE_ERROR );
00785         return;
00786     }
00787     // get some element info
00788     // We need to know how many of each element type we have.  We
00789     // are going to create an element block for each element type
00790     // that exists the mesh.  Block 1 will be tri3; block 2 will be
00791     // shell; block 3 will be tetra, and block 4 will be hex.
00792     const unsigned MAX_NODES = 27;
00793     const unsigned MIN_NODES = 2;
00794     int counts[MIXED][MAX_NODES + 1];
00795     memset( counts, 0, sizeof( counts ) );
00796 
00797     for( i = 0; i < myMesh->max_element_index(); ++i )
00798     {
00799         if( !myMesh->is_element_valid( i ) ) continue;
00800 
00801         EntityTopology type = myMesh->element_topology( i, err );MSQ_ERRRTN( err );
00802         unsigned nodes = myMesh->element_connectivity( i, err ).size();MSQ_ERRRTN( err );
00803         if( (unsigned)type >= MIXED || nodes < MIN_NODES || nodes > MAX_NODES )
00804         {
00805             MSQ_SETERR( err )( "Invalid element typology.", MsqError::INTERNAL_ERROR );
00806             return;
00807         }
00808         ++counts[type][nodes];
00809     }
00810 
00811     // Count number of required element blocks
00812     int block_count = 0;
00813     for( i = 0; i < MIXED; ++i )
00814         for( j = MIN_NODES; j < MAX_NODES; ++j )
00815             if( counts[i][j] ) ++block_count;
00816 
00817     // figure out if we have fixed nodes, if so, we need a nodeset
00818     int num_fixed_nodes = 0;
00819     for( i = 0; i < myMesh->max_vertex_index(); ++i )
00820     {
00821         bool fixed = myMesh->is_vertex_valid( i ) && myMesh->vertex_is_fixed( i, err );MSQ_ERRRTN( err );
00822         num_fixed_nodes += fixed;
00823     }
00824 
00825     // write doubles instead of floats
00826     int app_float_size  = sizeof( double );
00827     int file_float_size = sizeof( double );
00828     int exo_err         = 0;
00829 
00830     // Create the file.  If it exists, clobber it.  This could be dangerous.
00831     int file_id = ex_create( out_filename, EX_CLOBBER, &app_float_size, &file_float_size );
00832 
00833     // Make sure we opened the file correctly
00834     if( file_id < 0 )
00835     {
00836         MSQ_SETERR( err )( MsqError::FILE_ACCESS );
00837         return;
00838     }
00839 
00840     char title[MAX_LINE_LENGTH] = "Mesquite Export";
00841 
00842     size_t vert_count = myMesh->num_vertices();
00843     size_t elem_count = myMesh->num_elements();
00844 
00845     int ns_count = 0;
00846     if( num_fixed_nodes > 0 ) ns_count = 1;
00847     int ss_count = 0;
00848 
00849     // put the initial info about the file
00850     exo_err = ex_put_init( file_id, title, numCoords, vert_count, elem_count, block_count, ns_count, ss_count );
00851     if( exo_err < 0 )
00852     {
00853         MSQ_SETERR( err )( "Unable to initialize file data.", MsqError::IO_ERROR );
00854         return;
00855     }
00856 
00857     // Gather vertex coordinate data and write to file.
00858     std::vector< double > coords( vert_count * 3 );
00859     std::vector< double >::iterator x, y, z;
00860     x = coords.begin();
00861     y = x + vert_count;
00862     z = y + vert_count;
00863     for( i = 0; i < myMesh->max_vertex_index(); ++i )
00864     {
00865         if( !myMesh->is_vertex_valid( i ) ) continue;
00866 
00867         if( z == coords.end() )
00868         {
00869             MSQ_SETERR( err )( "Array overflow", MsqError::INTERNAL_ERROR );
00870             return;
00871         }
00872 
00873         Vector3D coords = myMesh->get_vertex_coords( i, err );MSQ_ERRRTN( err );
00874         *x = coords.x();
00875         ++x;
00876         *y = coords.y();
00877         ++y;
00878         *z = coords.z();
00879         ++z;
00880     }
00881     if( z != coords.end() )
00882     {
00883         MSQ_SETERR( err )( "Counter at incorrect number.", MsqError::INTERNAL_ERROR );
00884         return;
00885     }
00886     // put the coords
00887     exo_err = ex_put_coord( file_id, arrptr( coords ), &coords[vert_count], &coords[2 * vert_count] );
00888     if( exo_err < 0 )
00889     {
00890         MSQ_SETERR( err )( "Unable to put vertex coordinates in file.", MsqError::IO_ERROR );
00891         return;
00892     }
00893 
00894     // put the names of the coordinates
00895     char* coord_names[] = { "x", "y", "z" };
00896     exo_err             = ex_put_coord_names( file_id, coord_names );
00897 
00898     // Create element-type arrays indexed by EntityTopology
00899     const char* tri_name  = "TRI";
00900     const char* quad_name = "SHELL";
00901     const char* tet_name  = "TETRA";
00902     const char* hex_name  = "HEX";
00903     const char* wdg_name  = "WEDGE";
00904     const char* pyr_name  = "PYRAMID";
00905     const char* exo_names[MIXED];
00906     memset( exo_names, 0, sizeof( exo_names ) );
00907     exo_names[TRIANGLE]      = tri_name;
00908     exo_names[QUADRILATERAL] = quad_name;
00909     exo_names[TETRAHEDRON]   = tet_name;
00910     exo_names[HEXAHEDRON]    = hex_name;
00911     exo_names[PRISM]         = wdg_name;
00912     exo_names[PYRAMID]       = pyr_name;
00913     unsigned min_nodes[MIXED];
00914     memset( min_nodes, 0, sizeof( min_nodes ) );
00915     min_nodes[TRIANGLE]      = 3;
00916     min_nodes[QUADRILATERAL] = 4;
00917     min_nodes[TETRAHEDRON]   = 4;
00918     min_nodes[HEXAHEDRON]    = 8;
00919     min_nodes[PRISM]         = 6;
00920     min_nodes[PYRAMID]       = 5;
00921 
00922     // For each element type (topology and num nodes)
00923     int block_id = 0;
00924     char name_buf[16];
00925     int num_atts = 0;
00926     std::vector< int > conn;
00927     for( i = 0; i < MIXED; ++i )
00928     {
00929         for( j = MIN_NODES; j < MAX_NODES; ++j )
00930         {
00931             // Have any of this topo & count combination?
00932             if( !counts[i][j] ) continue;
00933 
00934             // Is a supported ExodusII type?
00935             if( !exo_names[i] )
00936             {
00937                 MSQ_SETERR( err )
00938                 ( MsqError::INVALID_STATE, "Element topology %d not supported by ExodusII", (int)i );
00939                 return;
00940             }
00941 
00942             // Construct ExodusII element name from topo & num nodes
00943             if( j == min_nodes[i] )
00944                 strcpy( name_buf, exo_names[i] );
00945             else
00946                 sprintf( name_buf, "%s%d", exo_names[i], (int)j );
00947 
00948             // Create element block
00949             ++block_id;
00950             exo_err = ex_put_elem_block( file_id, block_id, name_buf, counts[i][j], j, num_atts );
00951             if( exo_err < 0 )
00952             {
00953                 MSQ_SETERR( err )( "Error creating the tri block.", MsqError::IO_ERROR );
00954                 return;
00955             }
00956 
00957             // For each element
00958             conn.resize( counts[i][j] * j );
00959             std::vector< int >::iterator iter = conn.begin();
00960             for( k = 0; k < myMesh->max_element_index(); ++k )
00961             {
00962                 // If not correct topo, skip it.
00963                 if( !myMesh->is_element_valid( k ) || (unsigned)( myMesh->element_topology( k, err ) ) != i ) continue;MSQ_ERRRTN( err );
00964 
00965                 // If not correct number nodes, skip it
00966                 const std::vector< size_t >& elem_conn = myMesh->element_connectivity( k, err );MSQ_ERRRTN( err );
00967                 if( elem_conn.size() != j ) continue;
00968 
00969                 // Append element connectivity to list
00970                 for( std::vector< size_t >::const_iterator citer = elem_conn.begin(); citer != elem_conn.end();
00971                      ++citer, ++iter )
00972                 {
00973                     assert( iter != conn.end() );
00974                     *iter = *citer + 1;
00975                 }
00976             }
00977 
00978             // Make sure everything adds up
00979             if( iter != conn.end() )
00980             {
00981                 MSQ_SETERR( err )( MsqError::INTERNAL_ERROR );
00982                 return;
00983             }
00984 
00985             // Write element block connectivity
00986             exo_err = ex_put_elem_conn( file_id, block_id, arrptr( conn ) );
00987             if( exo_err < 0 )
00988             {
00989                 MSQ_SETERR( err )( "Error writing element connectivity.", MsqError::IO_ERROR );
00990                 return;
00991             }
00992         }
00993     }
00994 
00995     // Finally, mark boundary nodes
00996 
00997     if( num_fixed_nodes > 0 )
00998     {
00999         exo_err = ex_put_node_set_param( file_id, 111, num_fixed_nodes, 0 );
01000         if( exo_err < 0 )
01001         {
01002             MSQ_SETERR( err )( "Error while initializing node set.", MsqError::IO_ERROR );
01003             return;
01004         }
01005 
01006         int node_id = 0;
01007         std::vector< int > fixed_nodes( num_fixed_nodes );
01008         std::vector< int >::iterator iter = fixed_nodes.begin();
01009         for( i = 0; i < myMesh->max_vertex_index(); ++i )
01010         {
01011             if( !myMesh->is_vertex_valid( i ) ) continue;
01012             ++node_id;
01013 
01014             if( myMesh->vertex_is_fixed( i, err ) )
01015             {
01016                 if( iter == fixed_nodes.end() )
01017                 {
01018                     MSQ_SETERR( err )( MsqError::INTERNAL_ERROR );
01019                     return;
01020                 }
01021                 *iter = node_id;
01022                 ++iter;
01023             }
01024         }
01025 
01026         if( iter != fixed_nodes.end() )
01027         {
01028             MSQ_SETERR( err )( MsqError::INTERNAL_ERROR );
01029             return;
01030         }
01031 
01032         exo_err = ex_put_node_set( file_id, 111, arrptr( fixed_nodes ) );
01033         if( exo_err < 0 )
01034         {
01035             MSQ_SETERR( err )( "Error while writing node set.", MsqError::IO_ERROR );
01036             return;
01037         }
01038     }
01039     exo_err = ex_close( file_id );
01040     if( exo_err < 0 ) MSQ_SETERR( err )( "Error closing Exodus file.", MsqError::IO_ERROR );
01041 
01042 #endif
01043 }
01044 
01045 // Returns whether this mesh lies in a 2D or 3D coordinate system.
01046 
01047 int MeshImpl::get_geometric_dimension( MsqError& /*err*/ )
01048 {
01049     return numCoords;
01050 }
01051 
01052 void MeshImpl::get_all_elements( std::vector< ElementHandle >& elems, MsqError& err )
01053 {
01054     assert( sizeof( ElementHandle ) == sizeof( size_t ) );
01055     std::vector< size_t > temp;
01056     myMesh->all_elements( temp, err );MSQ_ERRRTN( err );
01057     elems.resize( temp.size() );
01058     if( !elems.empty() ) memcpy( arrptr( elems ), arrptr( temp ), sizeof( size_t ) * temp.size() );
01059 }
01060 
01061 void MeshImpl::get_all_vertices( std::vector< VertexHandle >& verts, MsqError& err )
01062 {
01063     assert( sizeof( VertexHandle ) == sizeof( size_t ) );
01064     std::vector< size_t > temp;
01065     myMesh->all_vertices( temp, err );MSQ_ERRRTN( err );
01066     verts.resize( temp.size() );
01067     if( !verts.empty() ) memcpy( arrptr( verts ), arrptr( temp ), sizeof( size_t ) * temp.size() );
01068 }
01069 
01070 // Returns a pointer to an iterator that iterates over the
01071 // set of all vertices in this mesh.  The calling code should
01072 // delete the returned iterator when it is finished with it.
01073 // If vertices are added or removed from the Mesh after obtaining
01074 // an iterator, the behavior of that iterator is undefined.
01075 VertexIterator* MeshImpl::vertex_iterator( MsqError& /*err*/ )
01076 {
01077     return new MeshImplVertIter( myMesh );
01078 }
01079 
01080 // Returns a pointer to an iterator that iterates over the
01081 // set of all top-level elements in this mesh.  The calling code should
01082 // delete the returned iterator when it is finished with it.
01083 // If elements are added or removed from the Mesh after obtaining
01084 // an iterator, the behavior of that iterator is undefined.
01085 ElementIterator* MeshImpl::element_iterator( MsqError& /*err*/ )
01086 {
01087     return new MeshImplElemIter( myMesh );
01088 }
01089 
01090 //************ Vertex Properties ********************
01091 // Returns true or false, indicating whether the vertex
01092 // is allowed to be repositioned.  True indicates that the vertex
01093 // is fixed and cannot be moved.  Note that this is a read-only
01094 // property; this flag can't be modified by users of the
01095 // Mesh interface.
01096 void MeshImpl::vertices_get_fixed_flag( const VertexHandle vert_array[], std::vector< bool >& on_bnd, size_t num_vtx,
01097                                         MsqError& err )
01098 {
01099     on_bnd.resize( num_vtx );
01100     for( size_t i = 0; i < num_vtx; ++i )
01101     {
01102         on_bnd[i] = myMesh->vertex_is_fixed( (size_t)vert_array[i], err );MSQ_ERRRTN( err );
01103     }
01104 }
01105 
01106 void MeshImpl::vertices_set_fixed_flag( const VertexHandle vert_array[], const std::vector< bool >& on_bnd,
01107                                         size_t num_vtx, MsqError& err )
01108 {
01109     assert( on_bnd.size() >= num_vtx );
01110     for( size_t i = 0; i < num_vtx; ++i )
01111     {
01112         myMesh->fix_vertex( (size_t)vert_array[i], on_bnd[i], err );MSQ_ERRRTN( err );
01113     }
01114 }
01115 void MeshImpl::vertices_get_slaved_flag( const VertexHandle vert_array[], std::vector< bool >& flags, size_t num_vtx,
01116                                          MsqError& err )
01117 {
01118     flags.resize( num_vtx );
01119     for( size_t i = 0; i < num_vtx; ++i )
01120     {
01121         flags[i] = myMesh->vertex_is_slaved( (size_t)vert_array[i], err );MSQ_ERRRTN( err );
01122     }
01123 }
01124 
01125 // Get/set location of a vertex
01126 void MeshImpl::vertices_get_coordinates( const Mesh::VertexHandle vert_array[], MsqVertex* coordinates, size_t num_vtx,
01127                                          MsqError& err )
01128 {
01129     for( size_t i = 0; i < num_vtx; ++i )
01130     {
01131         coordinates[i] = myMesh->get_vertex_coords( (size_t)vert_array[i], err );MSQ_ERRRTN( err );
01132     }
01133 }
01134 
01135 void MeshImpl::vertex_set_coordinates( VertexHandle vertex, const Vector3D& coordinates, MsqError& err )
01136 {
01137     myMesh->set_vertex_coords( (size_t)vertex, coordinates, err );MSQ_CHKERR( err );
01138 }
01139 
01140 // Each vertex has a byte-sized flag that can be used to store
01141 // flags.  This byte's value is neither set nor used by the mesh
01142 // implementation.  It is intended to be used by Mesquite algorithms.
01143 // Until a vertex's byte has been explicitly set, its value is 0.
01144 void MeshImpl::vertex_set_byte( VertexHandle vertex, unsigned char byte, MsqError& err )
01145 {
01146     vertices_set_byte( &vertex, &byte, 1, err );MSQ_CHKERR( err );
01147 }
01148 
01149 void MeshImpl::vertices_get_byte( const VertexHandle* vert_array, unsigned char* byte_array, size_t array_size,
01150                                   MsqError& err )
01151 {
01152     for( size_t i = 0; i < array_size; i++ )
01153     {
01154         byte_array[i] = myMesh->get_vertex_byte( (size_t)vert_array[i], err );MSQ_ERRRTN( err );
01155     }
01156 }
01157 
01158 // Retrieve the byte value for the specified vertex or vertices.
01159 // The byte value is 0 if it has not yet been set via one of the
01160 // *_set_byte() functions.
01161 void MeshImpl::vertex_get_byte( const VertexHandle vertex, unsigned char* byte, MsqError& err )
01162 {
01163     vertices_get_byte( &vertex, byte, 1, err );MSQ_CHKERR( err );
01164 }
01165 
01166 void MeshImpl::vertices_set_byte( const VertexHandle* vertex, const unsigned char* byte_array, size_t array_size,
01167                                   MsqError& err )
01168 {
01169     for( size_t i = 0; i < array_size; i++ )
01170     {
01171         myMesh->set_vertex_byte( (size_t)vertex[i], byte_array[i], err );MSQ_ERRRTN( err );
01172     }
01173 }
01174 
01175 template < typename T >
01176 struct cast_handle : public std::unary_function< size_t, T >
01177 {
01178     T operator()( size_t idx ) const
01179     {
01180         return reinterpret_cast< T >( idx );
01181     }
01182 };
01183 
01184 void MeshImpl::vertices_get_attached_elements( const VertexHandle* vertices, size_t num_vertices,
01185                                                std::vector< ElementHandle >& elements, std::vector< size_t >& offsets,
01186                                                MsqError& err )
01187 {
01188     elements.clear();
01189     offsets.clear();
01190     size_t prev_offset = 0;
01191     offsets.reserve( num_vertices + 1 );
01192     offsets.push_back( prev_offset );
01193     const VertexHandle* const vtx_end = vertices + num_vertices;
01194     for( ; vertices < vtx_end; ++vertices )
01195     {
01196         const std::vector< size_t >& adj = myMesh->vertex_adjacencies( (size_t)*vertices, err );MSQ_ERRRTN( err );
01197 
01198         prev_offset = prev_offset + adj.size();
01199         offsets.push_back( prev_offset );
01200 
01201         std::transform( adj.begin(), adj.end(), std::back_inserter( elements ), cast_handle< ElementHandle >() );
01202     }
01203 }
01204 
01205 void MeshImpl::elements_get_attached_vertices( const ElementHandle* elements, size_t num_elems,
01206                                                std::vector< VertexHandle >& vertices, std::vector< size_t >& offsets,
01207                                                MsqError& err )
01208 {
01209     vertices.clear();
01210     offsets.clear();
01211     size_t prev_offset = 0;
01212     offsets.reserve( num_elems + 1 );
01213     offsets.push_back( prev_offset );
01214     const ElementHandle* const elem_end = elements + num_elems;
01215     for( ; elements < elem_end; ++elements )
01216     {
01217         const std::vector< size_t >& conn = myMesh->element_connectivity( (size_t)*elements, err );MSQ_ERRRTN( err );
01218 
01219         prev_offset = prev_offset + conn.size();
01220         offsets.push_back( prev_offset );
01221 
01222         std::transform( conn.begin(), conn.end(), std::back_inserter( vertices ), cast_handle< VertexHandle >() );
01223     }
01224 }
01225 
01226 // Returns the topologies of the given entities.  The "entity_topologies"
01227 // array must be at least "num_elements" in size.
01228 void MeshImpl::elements_get_topologies( const ElementHandle* element_handle_array, EntityTopology* element_topologies,
01229                                         size_t num_elements, MsqError& err )
01230 {
01231     for( size_t i = 0; i < num_elements; i++ )
01232     {
01233         element_topologies[i] = myMesh->element_topology( (size_t)element_handle_array[i], err );MSQ_CHKERR( err );
01234     }
01235 }
01236 
01237 //**************** Memory Management ****************
01238 // Tells the mesh that the client is finished with a given
01239 // entity handle.
01240 void MeshImpl::release_entity_handles( const EntityHandle* /*handle_array*/, size_t /*num_handles*/, MsqError& /*err*/ )
01241 {
01242     // Do nothing
01243 }
01244 
01245 // Instead of deleting a Mesh when you think you are done,
01246 // call release().  In simple cases, the implementation could
01247 // just call the destructor.  More sophisticated implementations
01248 // may want to keep the Mesh object to live longer than Mesquite
01249 // is using it.
01250 void MeshImpl::release()
01251 {
01252     // delete this;
01253 }
01254 
01255 const char* const vtk_type_names[] = { "bit",          "char", "unsigned_char", "short", "unsigned_short", "int",
01256                                        "unsigned_int", "long", "unsigned_long", "float", "double",         0 };
01257 
01258 void MeshImpl::read_vtk( const char* filename, MsqError& err )
01259 {
01260     int major, minor;
01261     char vendor_string[257];
01262     size_t i;
01263 
01264     FILE* file = fopen( filename, "r" );
01265     if( !file )
01266     {
01267         MSQ_SETERR( err )( MsqError::FILE_ACCESS );
01268         return;
01269     }
01270 
01271     // Read file header
01272 
01273     if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
01274     {
01275         MSQ_SETERR( err )( MsqError::IO_ERROR );
01276         fclose( file );
01277         return;
01278     }
01279 
01280     if( !strchr( vendor_string, '\n' ) || 2 != sscanf( vendor_string, "# vtk DataFile Version %d.%d", &major, &minor ) )
01281     {
01282         MSQ_SETERR( err )( MsqError::FILE_FORMAT );
01283         fclose( file );
01284         return;
01285     }
01286 
01287     if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
01288     {
01289         MSQ_SETERR( err )( MsqError::IO_ERROR );
01290         fclose( file );
01291         return;
01292     }
01293 
01294     // VTK spec says this should not exceed 256 chars.
01295     if( !strchr( vendor_string, '\n' ) )
01296     {
01297         MSQ_SETERR( err )
01298         ( "Vendor string (line 2) exceeds 256 characters.", MsqError::PARSE_ERROR );
01299         fclose( file );
01300         return;
01301     }
01302 
01303     // Check file type
01304 
01305     FileTokenizer tokens( file );
01306     const char* const file_type_names[] = { "ASCII", "BINARY", 0 };
01307     int filetype                        = tokens.match_token( file_type_names, err );MSQ_ERRRTN( err );
01308     if( 2 == filetype )
01309     {
01310         MSQ_SETERR( err )
01311         ( "Cannot read BINARY VTK files -- use ASCII.", MsqError::NOT_IMPLEMENTED );
01312         return;
01313     }
01314 
01315     // Clear any existing data
01316     this->clear();
01317 
01318     const char* outer_block_names[] = { "DATASET", "FIELD", 0 };
01319     // Read the mesh
01320     // VTK docs are inconsistant with regard to whether or not
01321     // a field block should be specified as "FIELD" or "DATASET FIELD",
01322     // so allow both.
01323     while( !tokens.eof() )
01324     {
01325         int blocktype = tokens.match_token( outer_block_names, err );
01326         if( MSQ_CHKERR( err ) )
01327         {
01328             if( tokens.eof() )
01329             {
01330                 err.clear();
01331                 break;
01332             }
01333             else
01334                 return;
01335         }
01336 
01337         if( blocktype == 1 )
01338         {
01339             vtk_read_dataset( tokens, err );MSQ_ERRRTN( err );
01340         }
01341         else
01342         {
01343             vtk_read_field( tokens, err );MSQ_ERRRTN( err );
01344         }
01345     }
01346 
01347     // Make sure file actually contained some mesh
01348     if( myMesh->num_elements() == 0 )
01349     {
01350         MSQ_SETERR( err )( "File contained no mesh.", MsqError::PARSE_ERROR );
01351         return;
01352     }
01353 
01354     // There is no option for a 2-D mesh in VTK files.  Always 3
01355     numCoords = 3;
01356 
01357     // Convert tag data for fixed nodes to internal bitmap
01358     std::vector< bool > flags;
01359     tag_to_bool( "fixed", flags, err );MSQ_ERRRTN( err );
01360     if( !flags.empty() )
01361     {
01362         for( i = 0; i < myMesh->max_vertex_index(); ++i )
01363             myMesh->fix_vertex( i, flags[i], err );MSQ_ERRRTN( err );
01364     }
01365 
01366     flags.clear();
01367     tag_to_bool( "slaved", flags, err );MSQ_ERRRTN( err );
01368     if( !flags.empty() )
01369     {
01370         for( i = 0; i < myMesh->max_vertex_index(); ++i )
01371             myMesh->slave_vertex( i, flags[i], err );MSQ_ERRRTN( err );
01372     }
01373 }
01374 
01375 void MeshImpl::tag_to_bool( const char* tag_name, std::vector< bool >& values, MsqError& err )
01376 {
01377     // Convert tag data for fixed nodes to internal bitmap
01378     TagHandle handle = tag_get( tag_name, err );
01379     if( !handle || MSQ_CHKERR( err ) )
01380     {
01381         err.clear();
01382         values.clear();
01383         return;
01384     }
01385 
01386     size_t i;
01387     values.resize( myMesh->max_vertex_index(), false );
01388     const TagDescription& tag_desc = myTags->properties( (size_t)handle, err );MSQ_ERRRTN( err );
01389     bool havedata = myTags->tag_has_vertex_data( (size_t)handle, err );MSQ_ERRRTN( err );
01390     if( !havedata )
01391     {
01392         MSQ_SETERR( err )
01393         ( MsqError::FILE_FORMAT, "'%s' attribute on elements, not vertices", tag_name );
01394         return;
01395     }
01396 
01397     switch( tag_desc.type )
01398     {
01399         case BYTE: {
01400             char data;
01401             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01402             {
01403                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01404                 values[i] = !!data;
01405             }
01406             break;
01407         }
01408         case BOOL: {
01409             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01410             {
01411                 bool data;
01412                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01413                 values[i] = data;
01414             }
01415             break;
01416         }
01417         case INT: {
01418             int data;
01419             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01420             {
01421                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01422                 values[i] = !!data;
01423             }
01424             break;
01425         }
01426         case DOUBLE: {
01427             double data;
01428             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01429             {
01430                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01431                 values[i] = !!data;
01432             }
01433             break;
01434         }
01435         case HANDLE: {
01436             unsigned long data;
01437             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01438             {
01439                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01440                 values[i] = !!data;
01441             }
01442             break;
01443         }
01444         default:
01445             MSQ_SETERR( err )( MsqError::PARSE_ERROR, "'%s' attribute has invalid type", tag_name );
01446             return;
01447     }
01448 
01449     tag_delete( handle, err );
01450 }
01451 
01452 void MeshImpl::vtk_read_dataset( FileTokenizer& tokens, MsqError& err )
01453 {
01454     const char* const data_type_names[] = {
01455         "STRUCTURED_POINTS", "STRUCTURED_GRID", "UNSTRUCTURED_GRID", "POLYDATA", "RECTILINEAR_GRID", "FIELD", 0
01456     };
01457     int datatype = tokens.match_token( data_type_names, err );MSQ_ERRRTN( err );
01458 
01459     // Ignore FIELD data at beginning of DATASET. As far as I (J.Kraftcheck)
01460     // understand the VTK documentation, there should never be a FIELD block
01461     // inside a dataset, except in attribute data.  However, some app somewhere
01462     // writes them this way, so try to handle them.
01463     for( ;; )
01464     {
01465         tokens.match_token( "FIELD", err );
01466         if( err )
01467         {
01468             tokens.unget_token();
01469             err.clear();
01470             break;
01471         }
01472         vtk_read_field( tokens, err );MSQ_ERRRTN( err );
01473     }
01474 
01475     switch( datatype )
01476     {
01477         case 1:
01478             vtk_read_structured_points( tokens, err );
01479             break;
01480         case 2:
01481             vtk_read_structured_grid( tokens, err );
01482             break;
01483         case 3:
01484             vtk_read_unstructured_grid( tokens, err );
01485             break;
01486         case 4:
01487             vtk_read_polydata( tokens, err );
01488             break;
01489         case 5:
01490             vtk_read_rectilinear_grid( tokens, err );
01491             break;
01492         case 6:
01493             vtk_read_field( tokens, err );
01494             break;
01495     }
01496 
01497     // Read attribute data
01498     const char* const block_type_names[] = { "POINT_DATA", "CELL_DATA", "DATASET", 0 };
01499     int blocktype                        = 0;
01500     while( !tokens.eof() )
01501     {
01502         // get POINT_DATA or CELL_DATA
01503         int new_block_type = tokens.match_token( block_type_names, err );
01504         if( new_block_type == 3 )  // done reading attributes
01505         {
01506             tokens.unget_token();
01507             return;
01508         }
01509 
01510         if( err )
01511         {
01512             if( tokens.eof() )
01513             {
01514                 err.clear();
01515                 break;
01516             }
01517             // If next token was neither POINT_DATA nor CELL_DATA,
01518             // then there's another attribute under the current one.
01519             if( blocktype )
01520             {
01521                 tokens.unget_token();
01522                 err.clear();
01523             }
01524             else
01525             {
01526                 MSQ_ERRRTN( err );
01527             }
01528         }
01529         else
01530         {
01531             blocktype = new_block_type;
01532             long count;
01533             tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01534 
01535             if( blocktype == 1 && (unsigned long)count != myMesh->num_vertices() )
01536             {
01537                 MSQ_SETERR( err )
01538                 ( MsqError::PARSE_ERROR,
01539                   "Count inconsistent with number of vertices"
01540                   " at line %d.",
01541                   tokens.line_number() );
01542                 return;
01543             }
01544             else if( blocktype == 2 && (unsigned long)count != myMesh->num_elements() )
01545             {
01546                 MSQ_SETERR( err )
01547                 ( MsqError::PARSE_ERROR,
01548                   "Count inconsistent with number of elements"
01549                   " at line %d.",
01550                   tokens.line_number() );
01551                 return;
01552             }
01553         }
01554 
01555         if( blocktype == 1 )
01556             vtk_read_point_data( tokens, err );
01557         else
01558             vtk_read_cell_data( tokens, err );MSQ_ERRRTN( err );
01559     }
01560 }
01561 
01562 void MeshImpl::vtk_read_structured_points( FileTokenizer& tokens, MsqError& err )
01563 {
01564     long i, j, k;
01565     long dims[3];
01566     double origin[3], space[3];
01567 
01568     tokens.match_token( "DIMENSIONS", err );MSQ_ERRRTN( err );
01569     tokens.get_long_ints( 3, dims, err );MSQ_ERRRTN( err );
01570     tokens.get_newline( err );MSQ_ERRRTN( err );
01571 
01572     if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
01573     {
01574         MSQ_SETERR( err )
01575         ( MsqError::PARSE_ERROR, "Invalid dimension at line %d", tokens.line_number() );
01576         return;
01577     }
01578 
01579     tokens.match_token( "ORIGIN", err );MSQ_ERRRTN( err );
01580     tokens.get_doubles( 3, origin, err );MSQ_ERRRTN( err );
01581     tokens.get_newline( err );MSQ_ERRRTN( err );
01582 
01583     const char* const spacing_names[] = { "SPACING", "ASPECT_RATIO", 0 };
01584     tokens.match_token( spacing_names, err );MSQ_ERRRTN( err );
01585     tokens.get_doubles( 3, space, err );MSQ_ERRRTN( err );
01586     tokens.get_newline( err );MSQ_ERRRTN( err );
01587 
01588     myMesh->allocate_vertices( dims[0] * dims[1] * dims[2], err );MSQ_ERRRTN( err );
01589     size_t vtx = 0;
01590     Vector3D off( origin[0], origin[1], origin[2] );
01591     for( k = 0; k < dims[2]; ++k )
01592         for( j = 0; j < dims[1]; ++j )
01593             for( i = 0; i < dims[0]; ++i )
01594             {
01595                 myMesh->reset_vertex( vtx++, off + Vector3D( i * space[0], j * space[1], k * space[2] ), false, err );MSQ_ERRRTN( err );
01596             }
01597 
01598     vtk_create_structured_elems( dims, err );MSQ_ERRRTN( err );
01599 }
01600 
01601 void MeshImpl::vtk_read_structured_grid( FileTokenizer& tokens, MsqError& err )
01602 {
01603     long num_verts, dims[3];
01604 
01605     tokens.match_token( "DIMENSIONS", err );MSQ_ERRRTN( err );
01606     tokens.get_long_ints( 3, dims, err );MSQ_ERRRTN( err );
01607     tokens.get_newline( err );MSQ_ERRRTN( err );
01608 
01609     if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
01610     {
01611         MSQ_SETERR( err )
01612         ( MsqError::PARSE_ERROR, "Invalid dimension at line %d", tokens.line_number() );
01613         return;
01614     }
01615 
01616     tokens.match_token( "POINTS", err );MSQ_ERRRTN( err );
01617     tokens.get_long_ints( 1, &num_verts, err );MSQ_ERRRTN( err );
01618     tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01619     tokens.get_newline( err );MSQ_ERRRTN( err );
01620 
01621     if( num_verts != ( dims[0] * dims[1] * dims[2] ) )
01622     {
01623         MSQ_SETERR( err )
01624         ( MsqError::PARSE_ERROR,
01625           "Point count not consistent with dimensions "
01626           "at line %d",
01627           tokens.line_number() );
01628         return;
01629     }
01630 
01631     myMesh->allocate_vertices( num_verts, err );MSQ_ERRRTN( err );
01632     for( size_t vtx = 0; vtx < (size_t)num_verts; ++vtx )
01633     {
01634         Vector3D pos;
01635         tokens.get_doubles( 3, const_cast< double* >( pos.to_array() ), err );MSQ_ERRRTN( err );
01636         myMesh->reset_vertex( vtx, pos, false, err );MSQ_ERRRTN( err );
01637     }
01638 
01639     vtk_create_structured_elems( dims, err );MSQ_ERRRTN( err );
01640 }
01641 
01642 void MeshImpl::vtk_read_rectilinear_grid( FileTokenizer& tokens, MsqError& err )
01643 {
01644     int i, j, k;
01645     long dims[3];
01646     const char* labels[] = { "X_COORDINATES", "Y_COORDINATES", "Z_COORDINATES" };
01647     vector< double > coords[3];
01648 
01649     tokens.match_token( "DIMENSIONS", err );MSQ_ERRRTN( err );
01650     tokens.get_long_ints( 3, dims, err );MSQ_ERRRTN( err );
01651     tokens.get_newline( err );MSQ_ERRRTN( err );
01652 
01653     if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
01654     {
01655         MSQ_SETERR( err )
01656         ( MsqError::PARSE_ERROR, "Invalid dimension at line %d", tokens.line_number() );
01657         return;
01658     }
01659 
01660     for( i = 0; i < 3; i++ )
01661     {
01662         long count;
01663         tokens.match_token( labels[i], err );MSQ_ERRRTN( err );
01664         tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01665         tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01666         tokens.get_newline( err );MSQ_ERRRTN( err );
01667 
01668         if( count != dims[i] )
01669         {
01670             MSQ_SETERR( err )
01671             ( MsqError::PARSE_ERROR,
01672               "Coordinate count inconsistent with dimensions"
01673               " at line %d",
01674               tokens.line_number() );
01675             return;
01676         }
01677 
01678         coords[i].resize( count );
01679         tokens.get_doubles( count, &coords[i][0], err );MSQ_ERRRTN( err );
01680     }
01681 
01682     myMesh->allocate_vertices( dims[0] * dims[1] * dims[2], err );MSQ_ERRRTN( err );
01683     size_t vtx = 0;
01684     for( k = 0; k < dims[2]; ++k )
01685         for( j = 0; j < dims[1]; ++j )
01686             for( i = 0; i < dims[0]; ++i )
01687             {
01688                 myMesh->reset_vertex( vtx++, Vector3D( coords[0][i], coords[1][j], coords[2][k] ), false, err );MSQ_ERRRTN( err );
01689             }
01690 
01691     vtk_create_structured_elems( dims, err );MSQ_ERRRTN( err );
01692 }
01693 
01694 void MeshImpl::vtk_read_polydata( FileTokenizer& tokens, MsqError& err )
01695 {
01696     long num_verts;
01697     vector< int > connectivity;
01698     const char* const poly_data_names[] = { "VERTICES", "LINES", "POLYGONS", "TRIANGLE_STRIPS", 0 };
01699 
01700     tokens.match_token( "POINTS", err );MSQ_ERRRTN( err );
01701     tokens.get_long_ints( 1, &num_verts, err );MSQ_ERRRTN( err );
01702     tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01703     tokens.get_newline( err );MSQ_ERRRTN( err );
01704 
01705     if( num_verts < 1 )
01706     {
01707         MSQ_SETERR( err )
01708         ( MsqError::PARSE_ERROR, "Invalid point count at line %d", tokens.line_number() );
01709         return;
01710     }
01711 
01712     myMesh->allocate_vertices( num_verts, err );MSQ_ERRRTN( err );
01713     for( size_t vtx = 0; vtx < (size_t)num_verts; ++vtx )
01714     {
01715         Vector3D pos;
01716         tokens.get_doubles( 3, const_cast< double* >( pos.to_array() ), err );MSQ_ERRRTN( err );
01717         myMesh->reset_vertex( vtx, pos, false, err );MSQ_ERRRTN( err );
01718     }
01719 
01720     int poly_type = tokens.match_token( poly_data_names, err );MSQ_ERRRTN( err );
01721     switch( poly_type )
01722     {
01723         case 3:
01724             vtk_read_polygons( tokens, err );MSQ_ERRRTN( err );
01725             break;
01726         case 4:
01727             MSQ_SETERR( err )
01728             ( MsqError::NOT_IMPLEMENTED, "Unsupported type: triangle strips at line %d", tokens.line_number() );
01729             return;
01730         case 1:
01731         case 2:
01732             MSQ_SETERR( err )
01733             ( MsqError::NOT_IMPLEMENTED, "Entities of dimension < 2 at line %d", tokens.line_number() );
01734             return;
01735     }
01736 }
01737 
01738 void MeshImpl::vtk_read_polygons( FileTokenizer& tokens, MsqError& err )
01739 {
01740     long size[2];
01741 
01742     tokens.get_long_ints( 2, size, err );MSQ_ERRRTN( err );
01743     tokens.get_newline( err );MSQ_ERRRTN( err );
01744     myMesh->allocate_elements( size[0], err );MSQ_ERRRTN( err );
01745     std::vector< long > conn;
01746 
01747     for( int i = 0; i < size[0]; ++i )
01748     {
01749         long count;
01750         tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01751         conn.resize( count );
01752         tokens.get_long_ints( count, arrptr( conn ), err );MSQ_ERRRTN( err );
01753         myMesh->reset_element( i, conn, POLYGON, err );MSQ_ERRRTN( err );
01754     }
01755 }
01756 
01757 void MeshImpl::vtk_read_unstructured_grid( FileTokenizer& tokens, MsqError& err )
01758 {
01759     long i, num_verts, num_elems[2];
01760 
01761     tokens.match_token( "POINTS", err );MSQ_ERRRTN( err );
01762     tokens.get_long_ints( 1, &num_verts, err );MSQ_ERRRTN( err );
01763     tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01764     tokens.get_newline( err );MSQ_ERRRTN( err );
01765 
01766     if( num_verts < 1 )
01767     {
01768         MSQ_SETERR( err )
01769         ( MsqError::PARSE_ERROR, "Invalid point count at line %d", tokens.line_number() );
01770         return;
01771     }
01772 
01773     myMesh->allocate_vertices( num_verts, err );MSQ_ERRRTN( err );
01774     for( size_t vtx = 0; vtx < (size_t)num_verts; ++vtx )
01775     {
01776         Vector3D pos;
01777         tokens.get_doubles( 3, const_cast< double* >( pos.to_array() ), err );MSQ_ERRRTN( err );
01778         myMesh->reset_vertex( vtx, pos, false, err );MSQ_ERRRTN( err );
01779     }
01780 
01781     tokens.match_token( "CELLS", err );MSQ_ERRRTN( err );
01782     tokens.get_long_ints( 2, num_elems, err );MSQ_ERRRTN( err );
01783     tokens.get_newline( err );MSQ_ERRRTN( err );
01784 
01785     myMesh->allocate_elements( num_elems[0], err );MSQ_ERRRTN( err );
01786     std::vector< long > conn;
01787     for( i = 0; i < num_elems[0]; ++i )
01788     {
01789         long count;
01790         tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01791         conn.resize( count );
01792         tokens.get_long_ints( count, arrptr( conn ), err );MSQ_ERRRTN( err );
01793         myMesh->reset_element( i, conn, MIXED, err );MSQ_ERRRTN( err );
01794     }
01795 
01796     tokens.match_token( "CELL_TYPES", err );MSQ_ERRRTN( err );
01797     tokens.get_long_ints( 1, &num_elems[1], err );MSQ_ERRRTN( err );
01798     tokens.get_newline( err );MSQ_ERRRTN( err );
01799 
01800     if( num_elems[0] != num_elems[1] )
01801     {
01802         MSQ_SETERR( err )
01803         ( MsqError::PARSE_ERROR,
01804           "Number of element types does not match number of elements"
01805           "at line %d",
01806           tokens.line_number() );
01807         return;
01808     }
01809 
01810     std::vector< size_t > tconn;
01811     for( i = 0; i < num_elems[0]; ++i )
01812     {
01813         long type;
01814         size_t size;
01815         tokens.get_long_ints( 1, &type, err );MSQ_ERRRTN( err );
01816 
01817         // Check if type is a valid value
01818         const VtkTypeInfo* info = VtkTypeInfo::find_type( type, err );
01819         if( err || !info || ( !info->numNodes && type != POLYGON ) )
01820         {
01821             MSQ_SETERR( err )
01822             ( MsqError::PARSE_ERROR, "Invalid cell type %ld at line %d.", type, tokens.line_number() );
01823             return;
01824         }
01825         // Check if Mesquite supports the type
01826         if( info->msqType == MIXED )
01827         {
01828             MSQ_SETERR( err )
01829             ( MsqError::UNSUPPORTED_ELEMENT, "Unsupported cell type %ld (%s) at line %d.", type, info->name,
01830               tokens.line_number() );
01831             return;
01832         }
01833 
01834         // If node-ordering is not the same as exodus...
01835         if( info->vtkConnOrder )
01836         {
01837             size = myMesh->element_connectivity( i, err ).size();MSQ_ERRRTN( err );
01838             if( info->numNodes != size )
01839             {
01840                 MSQ_SETERR( err )
01841                 ( MsqError::UNSUPPORTED_ELEMENT, "Cell type %ld (%s) for element with %d nodes at Line %d", type,
01842                   info->name, (int)size, tokens.line_number() );
01843                 return;
01844             }
01845 
01846             tconn.resize( size );
01847             const std::vector< size_t >& pconn = myMesh->element_connectivity( i, err );MSQ_ERRRTN( err );
01848             for( size_t j = 0; j < size; ++j )
01849             {
01850                 tconn[j] = pconn[info->vtkConnOrder[j]];
01851             }
01852 
01853             myMesh->reset_element( i, tconn, info->msqType, err );MSQ_ERRRTN( err );
01854         }
01855         // Othewise (if node ordering is the same), just set the type.
01856         else
01857         {
01858             myMesh->element_topology( i, info->msqType, err );MSQ_ERRRTN( err );
01859         }
01860     }  // for(i)
01861 }
01862 
01863 void MeshImpl::vtk_create_structured_elems( const long* dims, MsqError& err )
01864 {
01865     // NOTE: this should be work fine for edges also if
01866     //      Mesquite ever supports them.  Just add the
01867     //      type for dimension 1 to the switch statement.
01868 
01869     // int non_zero[3] = {0,0,0};  // True if dim > 0 for x, y, z respectively
01870     long elem_dim  = 0;           // Element dimension (2->quad, 3->hex)
01871     long num_elems = 1;           // Total number of elements
01872     long vert_per_elem;           // Element connectivity length
01873     long edims[3] = { 1, 1, 1 };  // Number of elements in each grid direction
01874 
01875     // Populate above data
01876     for( int d = 0; d < 3; d++ )
01877         if( dims[d] > 1 )
01878         {
01879             // non_zero[elem_dim] = d;
01880             edims[d] = dims[d] - 1;
01881             num_elems *= edims[d];
01882             elem_dim++;
01883         }
01884     vert_per_elem = 1 << elem_dim;
01885 
01886     // Get element type from element dimension
01887     EntityTopology type;
01888     switch( elem_dim )
01889     {
01890             // case 1: type = EDGE;          break;
01891         case 2:
01892             type = QUADRILATERAL;
01893             break;
01894         case 3:
01895             type = HEXAHEDRON;
01896             break;
01897         default:
01898             MSQ_SETERR( err )
01899             ( "Cannot create structured mesh with elements "
01900               "of dimension < 2 or > 3.",
01901               MsqError::NOT_IMPLEMENTED );
01902             return;
01903     }
01904 
01905     // Allocate storage for elements
01906     myMesh->allocate_elements( num_elems, err );MSQ_ERRRTN( err );
01907 
01908     // Offsets of element vertices in grid relative to corner closest to origin
01909     long k                = dims[0] * dims[1];
01910     const long corners[8] = { 0, 1, 1 + dims[0], dims[0], k, k + 1, k + 1 + dims[0], k + dims[0] };
01911 
01912     // Populate element list
01913     std::vector< size_t > conn( vert_per_elem );
01914     size_t elem_idx = 0;
01915     for( long z = 0; z < edims[2]; ++z )
01916         for( long y = 0; y < edims[1]; ++y )
01917             for( long x = 0; x < edims[0]; ++x )
01918             {
01919                 const long index = x + y * dims[0] + z * ( dims[0] * dims[1] );
01920                 for( long j = 0; j < vert_per_elem; ++j )
01921                     conn[j] = index + corners[j];
01922                 myMesh->reset_element( elem_idx++, conn, type, err );MSQ_ERRRTN( err );
01923             }
01924 }
01925 
01926 void* MeshImpl::vtk_read_field_data( FileTokenizer& tokens, size_t count, size_t num_fields,
01927                                      const std::string& field_name, TagDescription& tag, MsqError& err )
01928 {
01929     tag.member = tokens.get_string( err );
01930     MSQ_ERRZERO( err );
01931     // If field is a struct containing multiple members, make
01932     // tag name the concatentation of the field name and member
01933     // name.
01934     if( num_fields > 1 )
01935     {
01936         tag.name = field_name;
01937         tag.name += " ";
01938         tag.name += tag.member;
01939         tag.member.clear();
01940     }
01941     // If field contains only one member, then make the tag
01942     // name be the field name, and store the member name to
01943     // preserve it for subsequent writes.
01944     else
01945     {
01946         tag.name = field_name;
01947     }
01948 
01949     // Get tuple size and count from the file.
01950     long sizes[2];
01951     tokens.get_long_ints( 2, sizes, err );
01952     MSQ_ERRZERO( err );
01953     if( sizes[0] < 1 )
01954     {
01955         MSQ_SETERR( err )
01956         ( MsqError::PARSE_ERROR, "Invalid tuple size (%ld) for field data %s at line %d\n", sizes[0], tag.name.c_str(),
01957           tokens.line_number() );
01958         return 0;
01959     }
01960     if( sizes[1] < 0 || ( count && (size_t)sizes[1] != count ) )
01961     {
01962         MSQ_SETERR( err )
01963         ( MsqError::PARSE_ERROR,
01964           "Invalid field data length at line %d.  "
01965           "Cannot map %lu tuples to  %ld entities.\n",
01966           tokens.line_number(), (unsigned long)count, sizes[1] );
01967         return 0;
01968     }
01969 
01970     int type = tokens.match_token( vtk_type_names, err );
01971     MSQ_ERRZERO( err );
01972     void* result = vtk_read_typed_data( tokens, type, sizes[0], sizes[1], tag, err );MSQ_CHKERR( err );
01973     return result;
01974 }
01975 
01976 void MeshImpl::vtk_read_field( FileTokenizer& tokens, MsqError& err )
01977 {
01978     // This is not supported yet.
01979     // Parse the data but throw it away because
01980     // Mesquite has no internal representation for it.
01981 
01982     std::string name = tokens.get_string( err );MSQ_ERRRTN( err );
01983     int count;
01984     tokens.get_integers( 1, &count, err );MSQ_ERRRTN( err );
01985 
01986     for( int i = 0; i < count; i++ )
01987     {
01988         TagDescription tag;
01989         void* ptr = vtk_read_field_data( tokens, 0, 1, "", tag, err );
01990         if( ptr ) free( ptr );MSQ_ERRRTN( err );
01991     }
01992 }
01993 
01994 void* MeshImpl::vtk_read_attrib_data( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
01995 {
01996     const char* const type_names[] = { "SCALARS", "COLOR_SCALARS", "VECTORS", "NORMALS", "TEXTURE_COORDINATES",
01997                                        "TENSORS", "FIELD",
01998                                        // Some buggy VTK files have empty CELL_DATA/POINT_DATA
01999                                        // blocks Try to allow them by checking for possible other
02000                                        // tokens indicating the end of the block
02001                                        "POINT_DATA", "CELL_DATA", "DATASET", 0 };
02002 
02003     int type = tokens.match_token( type_names, err );
02004     MSQ_ERRZERO( err );
02005     if( type > 7 )  // empty CELL_DATA/POINT_DATA block.
02006     {
02007         tag.vtkType = TagDescription::NONE;
02008         tokens.unget_token();
02009         return 0;
02010     }
02011 
02012     const char* name = tokens.get_string( err );
02013     MSQ_ERRZERO( err );
02014     tag.name = name;
02015 
02016     void* data = 0;
02017     switch( type )
02018     {
02019         case 1:
02020             data        = vtk_read_scalar_attrib( tokens, count, tag, err );
02021             tag.vtkType = TagDescription::SCALAR;
02022             break;
02023         case 2:
02024             data        = vtk_read_color_attrib( tokens, count, tag, err );
02025             tag.vtkType = TagDescription::COLOR;
02026             break;
02027         case 3:
02028             data        = vtk_read_vector_attrib( tokens, count, tag, err );
02029             tag.vtkType = TagDescription::VECTOR;
02030             break;
02031         case 4:
02032             data        = vtk_read_vector_attrib( tokens, count, tag, err );
02033             tag.vtkType = TagDescription::NORMAL;
02034             break;
02035         case 5:
02036             data        = vtk_read_texture_attrib( tokens, count, tag, err );
02037             tag.vtkType = TagDescription::TEXTURE;
02038             break;
02039         case 6:
02040             data        = vtk_read_tensor_attrib( tokens, count, tag, err );
02041             tag.vtkType = TagDescription::TENSOR;
02042             break;
02043         case 7:
02044             data        = 0;
02045             tag.vtkType = TagDescription::FIELD;
02046             break;
02047     }
02048 
02049     return data;
02050 }
02051 
02052 void MeshImpl::vtk_read_point_data( FileTokenizer& tokens, MsqError& err )
02053 {
02054     TagDescription tag;
02055     void* data = vtk_read_attrib_data( tokens, myMesh->num_vertices(), tag, err );
02056     if( data )  // normal attribute
02057     {
02058         vtk_store_point_data( data, tag, err );
02059         free( data );MSQ_CHKERR( err );
02060         return;
02061         ;
02062     }
02063     else if( tag.vtkType == TagDescription::FIELD )
02064     {
02065         std::string field = tag.name;
02066         int field_count;
02067         tokens.get_integers( 1, &field_count, err );MSQ_ERRRTN( err );
02068         for( int i = 0; i < field_count; ++i )
02069         {
02070             data = vtk_read_field_data( tokens, myMesh->num_vertices(), field_count, field, tag, err );MSQ_ERRRTN( err );
02071             vtk_store_point_data( data, tag, err );
02072             free( data );MSQ_ERRRTN( err );
02073         }
02074     }
02075 }
02076 
02077 void MeshImpl::vtk_store_point_data( const void* data, TagDescription& tag, MsqError& err )
02078 {
02079     size_t tag_handle = myTags->handle( tag.name, err );MSQ_ERRRTN( err );
02080     if( !tag_handle )
02081     {
02082         tag_handle = myTags->create( tag, 0, err );MSQ_ERRRTN( err );
02083     }
02084     else
02085     {
02086         const TagDescription& desc = myTags->properties( tag_handle, err );MSQ_ERRRTN( err );
02087         if( desc != tag )
02088         {
02089             MSQ_SETERR( err )
02090             ( MsqError::PARSE_ERROR,
02091               "Inconsistent types between element "
02092               "and vertex attributes named \"%s\"\n",
02093               tag.name.c_str() );
02094             return;
02095         }
02096     }
02097 
02098     std::vector< size_t > vertex_handles;
02099     myMesh->all_vertices( vertex_handles, err );MSQ_ERRRTN( err );
02100     if( !vertex_handles.empty() )
02101         myTags->set_vertex_data( tag_handle, vertex_handles.size(), arrptr( vertex_handles ), data, err );MSQ_ERRRTN( err );
02102 }
02103 
02104 void MeshImpl::vtk_read_cell_data( FileTokenizer& tokens, MsqError& err )
02105 {
02106     TagDescription tag;
02107     void* data = vtk_read_attrib_data( tokens, myMesh->num_elements(), tag, err );
02108     if( data )  // normal attribute
02109     {
02110         vtk_store_cell_data( data, tag, err );
02111         free( data );MSQ_CHKERR( err );
02112         return;
02113         ;
02114     }
02115     else if( tag.vtkType == TagDescription::FIELD )
02116     {
02117         std::string field = tag.name;
02118         int field_count;
02119         tokens.get_integers( 1, &field_count, err );MSQ_ERRRTN( err );
02120         for( int i = 0; i < field_count; ++i )
02121         {
02122             data = vtk_read_field_data( tokens, myMesh->num_elements(), field_count, field, tag, err );MSQ_ERRRTN( err );
02123             vtk_store_cell_data( data, tag, err );
02124             free( data );MSQ_ERRRTN( err );
02125         }
02126     }
02127 }
02128 
02129 void MeshImpl::vtk_store_cell_data( const void* data, TagDescription& tag, MsqError& err )
02130 {
02131     size_t tag_handle = myTags->handle( tag.name, err );MSQ_ERRRTN( err );
02132     if( !tag_handle )
02133     {
02134         tag_handle = myTags->create( tag, 0, err );MSQ_ERRRTN( err );
02135     }
02136     else
02137     {
02138         const TagDescription& desc = myTags->properties( tag_handle, err );MSQ_ERRRTN( err );
02139         if( desc != tag )
02140         {
02141             MSQ_SETERR( err )
02142             ( MsqError::PARSE_ERROR,
02143               "Inconsistent types between element "
02144               "and vertex attributes named \"%s\"\n",
02145               tag.name.c_str() );
02146             return;
02147         }
02148     }
02149 
02150     std::vector< size_t > element_handles;
02151     myMesh->all_elements( element_handles, err );MSQ_ERRRTN( err );
02152     if( !element_handles.empty() )
02153         myTags->set_element_data( tag_handle, element_handles.size(), arrptr( element_handles ), data, err );MSQ_ERRRTN( err );
02154 }
02155 
02156 void* MeshImpl::vtk_read_typed_data( FileTokenizer& tokens, int type, size_t per_elem, size_t num_elem,
02157                                      TagDescription& tag, MsqError& err )
02158 {
02159     void* data_ptr;
02160     size_t count = per_elem * num_elem;
02161     switch( type )
02162     {
02163         case 1:
02164             tag.size = per_elem * sizeof( bool );
02165             tag.type = BOOL;
02166             data_ptr = malloc( num_elem * tag.size );
02167             tokens.get_booleans( count, (bool*)data_ptr, err );
02168             break;
02169         case 2:
02170         case 3:
02171         case 4:
02172         case 5:
02173         case 6:
02174         case 7:
02175             tag.size = per_elem * sizeof( int );
02176             tag.type = INT;
02177             data_ptr = malloc( num_elem * tag.size );
02178             tokens.get_integers( count, (int*)data_ptr, err );
02179             break;
02180         case 8:
02181         case 9:
02182             // this is a bit of a hack since MeshImpl doesn't have a LONG type (HANDLE is used by
02183             // ParallelMesh for long)
02184             tag.size = per_elem * sizeof( size_t );
02185             assert( sizeof( long ) == sizeof( size_t ) );
02186             assert( sizeof( long ) == sizeof( void* ) );
02187             tag.type = HANDLE;
02188             data_ptr = malloc( num_elem * tag.size );
02189             tokens.get_long_ints( count, (long*)data_ptr, err );
02190             break;
02191         case 10:
02192         case 11:
02193             tag.size = per_elem * sizeof( double );
02194             tag.type = DOUBLE;
02195             data_ptr = malloc( num_elem * tag.size );
02196             tokens.get_doubles( count, (double*)data_ptr, err );
02197             break;
02198         default:
02199             MSQ_SETERR( err )( "Invalid data type", MsqError::INVALID_ARG );
02200             return 0;
02201     }
02202 
02203     if( MSQ_CHKERR( err ) )
02204     {
02205         free( data_ptr );
02206         return 0;
02207     }
02208 
02209     return data_ptr;
02210 }
02211 
02212 void* MeshImpl::vtk_read_scalar_attrib( FileTokenizer& tokens, long count, TagDescription& desc, MsqError& err )
02213 {
02214     if( !count ) return 0;
02215 
02216     int type = tokens.match_token( vtk_type_names, err );
02217     MSQ_ERRZERO( err );
02218 
02219     long size;
02220     const char* tok = tokens.get_string( err );
02221     MSQ_ERRZERO( err );
02222     const char* end = 0;
02223     size            = strtol( tok, (char**)&end, 0 );
02224     if( *end )
02225     {
02226         size = 1;
02227         tokens.unget_token();
02228     }
02229 
02230     // VTK spec says cannot be greater than 4--do we care?
02231     if( size < 1 || size > 4 )
02232     {
02233         MSQ_SETERR( err )
02234         ( MsqError::PARSE_ERROR,
02235           "Scalar count out of range [1,4]"
02236           " at line %d",
02237           tokens.line_number() );
02238         return 0;
02239     }
02240 
02241     tokens.match_token( "LOOKUP_TABLE", err );
02242     MSQ_ERRZERO( err );
02243     tok = tokens.get_string( err );
02244     MSQ_ERRZERO( err );
02245 
02246     // If no lookup table, just read and return the data
02247     if( !strcmp( tok, "default" ) )
02248     {
02249         void* ptr = vtk_read_typed_data( tokens, type, size, count, desc, err );
02250         MSQ_ERRZERO( err );
02251         return ptr;
02252     }
02253 
02254     // If we got this far, then the data has a lookup
02255     // table.  First read the lookup table and convert
02256     // to integers.
02257     string name = tok;
02258     vector< long > table( size * count );
02259     if( type > 0 && type < 10 )  // Is an integer type
02260     {
02261         tokens.get_long_ints( table.size(), arrptr( table ), err );
02262         MSQ_ERRZERO( err );
02263     }
02264     else  // Is a real-number type
02265     {
02266         for( std::vector< long >::iterator iter = table.begin(); iter != table.end(); ++iter )
02267         {
02268             double data;
02269             tokens.get_doubles( 1, &data, err );
02270             MSQ_ERRZERO( err );
02271 
02272             *iter = (long)data;
02273             if( (double)*iter != data )
02274             {
02275                 MSQ_SETERR( err )
02276                 ( MsqError::PARSE_ERROR, "Invalid lookup index (%.0f) at line %d", data, tokens.line_number() );
02277                 return 0;
02278             }
02279         }
02280     }
02281 
02282     // Now read the data - must be float RGBA color triples
02283 
02284     long table_count;
02285     tokens.match_token( "LOOKUP_TABLE", err );
02286     MSQ_ERRZERO( err );
02287     tokens.match_token( name.c_str(), err );
02288     MSQ_ERRZERO( err );
02289     tokens.get_long_ints( 1, &table_count, err );
02290     MSQ_ERRZERO( err );
02291 
02292     vector< float > table_data( table_count * 4 );
02293     tokens.get_floats( table_data.size(), arrptr( table_data ), err );
02294     MSQ_ERRZERO( err );
02295 
02296     // Create list from indexed data
02297 
02298     float* data      = (float*)malloc( sizeof( float ) * count * size * 4 );
02299     float* data_iter = data;
02300     for( std::vector< long >::iterator idx = table.begin(); idx != table.end(); ++idx )
02301     {
02302         if( *idx < 0 || *idx >= table_count )
02303         {
02304             MSQ_SETERR( err )
02305             ( MsqError::PARSE_ERROR, "LOOKUP_TABLE index %ld out of range.", *idx );
02306             free( data );
02307             return 0;
02308         }
02309 
02310         for( int i = 0; i < 4; i++ )
02311         {
02312             *data_iter = table_data[4 * *idx + i];
02313             ++data_iter;
02314         }
02315     }
02316 
02317     desc.size = size * 4 * sizeof( float );
02318     desc.type = DOUBLE;
02319     return data;
02320 }
02321 
02322 void* MeshImpl::vtk_read_color_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02323 {
02324     long size;
02325     tokens.get_long_ints( 1, &size, err );
02326     MSQ_ERRZERO( err );
02327 
02328     if( size < 1 )
02329     {
02330         MSQ_SETERR( err )
02331         ( MsqError::PARSE_ERROR, "Invalid size (%ld) at line %d", size, tokens.line_number() );
02332         return 0;
02333     }
02334 
02335     float* data = (float*)malloc( sizeof( float ) * count * size );
02336     tokens.get_floats( count * size, data, err );
02337     if( MSQ_CHKERR( err ) )
02338     {
02339         free( data );
02340         return 0;
02341     }
02342 
02343     tag.size = size * sizeof( float );
02344     tag.type = DOUBLE;
02345     return data;
02346 }
02347 
02348 void* MeshImpl::vtk_read_vector_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02349 {
02350     int type = tokens.match_token( vtk_type_names, err );
02351     MSQ_ERRZERO( err );
02352 
02353     void* result = vtk_read_typed_data( tokens, type, 3, count, tag, err );
02354     MSQ_ERRZERO( err );
02355     return result;
02356 }
02357 
02358 void* MeshImpl::vtk_read_texture_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02359 {
02360     int type, dim;
02361     tokens.get_integers( 1, &dim, err );
02362     MSQ_ERRZERO( err );
02363     type = tokens.match_token( vtk_type_names, err );
02364     MSQ_ERRZERO( err );
02365 
02366     if( dim < 1 || dim > 3 )
02367     {
02368         MSQ_SETERR( err )
02369         ( MsqError::PARSE_ERROR, "Invalid dimension (%d) at line %d.", dim, tokens.line_number() );
02370         return 0;
02371     }
02372 
02373     void* result = vtk_read_typed_data( tokens, type, dim, count, tag, err );
02374     MSQ_ERRZERO( err );
02375     return result;
02376 }
02377 
02378 void* MeshImpl::vtk_read_tensor_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02379 {
02380     int type = tokens.match_token( vtk_type_names, err );
02381     MSQ_ERRZERO( err );
02382 
02383     void* result = vtk_read_typed_data( tokens, type, 9, count, tag, err );
02384     MSQ_ERRZERO( err );
02385     return result;
02386 }
02387 
02388 void MeshImpl::vtk_write_attrib_data( std::ostream& file, const TagDescription& desc, const void* data, size_t count,
02389                                       MsqError& err ) const
02390 {
02391     // srkenno@sandia.gov: we now allow this type to be able to write e.g. GLOBAL_ID for parallel
02392     // meshes
02393     /*
02394     if (desc.type == HANDLE)
02395     {
02396       MSQ_SETERR(err)("Cannot write HANDLE tag data to VTK file.",
02397                       MsqError::FILE_FORMAT);
02398       return;
02399       }*/
02400 
02401     TagDescription::VtkType vtk_type = desc.vtkType;
02402     unsigned vlen                    = desc.size / MeshImplTags::size_from_tag_type( desc.type );
02403     // guess one from data length if not set
02404     if( vtk_type == TagDescription::NONE )
02405     {
02406         switch( vlen )
02407         {
02408             default:
02409                 vtk_type = TagDescription::SCALAR;
02410                 break;
02411             case 3:
02412                 vtk_type = TagDescription::VECTOR;
02413                 break;
02414             case 9:
02415                 vtk_type = TagDescription::TENSOR;
02416                 break;
02417                 return;
02418         }
02419     }
02420 
02421     // srkenno@sandia.gov: from class Mesh, the typenames below should correspond in order...
02422     // enum TagType { BYTE, BOOL, INT, DOUBLE, HANDLE };
02423 
02424     const char* const typenames[] = { "unsigned_char", "bit", "int", "double", "unsigned_long" };
02425     std::string field, member;
02426 
02427     int num_per_line;
02428     switch( vtk_type )
02429     {
02430         case TagDescription::SCALAR:
02431             num_per_line = vlen;
02432             file << "SCALARS " << desc.name << " " << typenames[desc.type] << " " << vlen << "\n";
02433             file << "LOOKUP_TABLE default\n";
02434             break;
02435         case TagDescription::COLOR:
02436             num_per_line = vlen;
02437             file << "COLOR_SCALARS " << desc.name << " " << vlen << "\n";
02438             break;
02439         case TagDescription::VECTOR:
02440             num_per_line = 3;
02441             if( vlen != 3 )
02442             {
02443                 MSQ_SETERR( err )
02444                 ( MsqError::INTERNAL_ERROR, "Tag \"%s\" is labeled as a VTK vector attribute but has %u values.",
02445                   desc.name.c_str(), vlen );
02446                 return;
02447             }
02448             file << "VECTORS " << desc.name << " " << typenames[desc.type] << "\n";
02449             break;
02450         case TagDescription::NORMAL:
02451             num_per_line = 3;
02452             if( vlen != 3 )
02453             {
02454                 MSQ_SETERR( err )
02455                 ( MsqError::INTERNAL_ERROR, "Tag \"%s\" is labeled as a VTK normal attribute but has %u values.",
02456                   desc.name.c_str(), vlen );
02457                 return;
02458             }
02459             file << "NORMALS " << desc.name << " " << typenames[desc.type] << "\n";
02460             break;
02461         case TagDescription::TEXTURE:
02462             num_per_line = vlen;
02463             file << "TEXTURE_COORDINATES " << desc.name << " " << typenames[desc.type] << " " << vlen << "\n";
02464             break;
02465         case TagDescription::TENSOR:
02466             num_per_line = 3;
02467             if( vlen != 9 )
02468             {
02469                 MSQ_SETERR( err )
02470                 ( MsqError::INTERNAL_ERROR, "Tag \"%s\" is labeled as a VTK tensor attribute but has %u values.",
02471                   desc.name.c_str(), vlen );
02472                 return;
02473             }
02474             file << "TENSORS " << desc.name << " " << typenames[desc.type] << "\n";
02475             break;
02476         case TagDescription::FIELD:
02477             num_per_line = vlen;
02478             get_field_names( desc, field, member, err );MSQ_ERRRTN( err );
02479             file << member << " " << vlen << " " << count << " " << typenames[desc.type] << "\n";
02480             break;
02481         default:
02482             MSQ_SETERR( err )( "Unknown VTK attribute type for tag.", MsqError::INTERNAL_ERROR );
02483             return;
02484     }
02485 
02486     size_t i, total = count * vlen;
02487     char* space = new char[num_per_line];
02488     memset( space, ' ', num_per_line );
02489     space[num_per_line - 1]    = '\n';
02490     const unsigned char* odata = (const unsigned char*)data;
02491     const bool* bdata          = (const bool*)data;
02492     const int* idata           = (const int*)data;
02493     const long* ldata          = (const long*)data;
02494     const double* ddata        = (const double*)data;
02495     switch( desc.type )
02496     {
02497         case BYTE:
02498             for( i = 0; i < total; ++i )
02499                 file << (unsigned int)odata[i] << space[i % num_per_line];
02500             break;
02501         case BOOL:
02502             for( i = 0; i < total; ++i )
02503                 file << ( bdata[i] ? '1' : '0' ) << space[i % num_per_line];
02504             break;
02505         case INT:
02506             for( i = 0; i < total; ++i )
02507                 file << idata[i] << space[i % num_per_line];
02508             break;
02509         case DOUBLE:
02510             for( i = 0; i < total; ++i )
02511                 file << ddata[i] << space[i % num_per_line];
02512             break;
02513         case HANDLE:
02514             for( i = 0; i < total; ++i )
02515                 file << ldata[i] << space[i % num_per_line];
02516             break;
02517         default:
02518             MSQ_SETERR( err )( "Unknown tag type.", MsqError::INTERNAL_ERROR );
02519     }
02520     delete[] space;
02521 }
02522 
02523 /**************************************************************************
02524  *                               TAGS
02525  **************************************************************************/
02526 
02527 TagHandle MeshImpl::tag_create( const std::string& name, TagType type, unsigned length, const void* defval,
02528                                 MsqError& err )
02529 {
02530     TagDescription::VtkType vtype;
02531     std::string field;
02532     switch( length )
02533     {
02534         case 1:
02535             vtype = TagDescription::SCALAR;
02536             break;
02537         case 3:
02538             vtype = TagDescription::VECTOR;
02539             break;
02540         case 9:
02541             vtype = TagDescription::TENSOR;
02542             break;
02543         default:
02544             vtype = TagDescription::FIELD;
02545             field = MESQUITE_FIELD_TAG;
02546             break;
02547     }
02548 
02549     // If tag name contains a space, assume the tag name
02550     // is a concatenation of the VTK field and member names.
02551     if( vtype != TagDescription::FIELD && name.find( " " ) != std::string::npos ) vtype = TagDescription::FIELD;
02552 
02553     size_t size = MeshImplTags::size_from_tag_type( type );
02554     TagDescription desc( name, type, vtype, length * size, field );
02555     size_t index = myTags->create( desc, defval, err );
02556     MSQ_ERRZERO( err );
02557     return (TagHandle)index;
02558 }
02559 
02560 void MeshImpl::tag_delete( TagHandle handle, MsqError& err )
02561 {
02562     myTags->destroy( (size_t)handle, err );MSQ_CHKERR( err );
02563 }
02564 
02565 TagHandle MeshImpl::tag_get( const std::string& name, MsqError& err )
02566 {
02567     size_t index = myTags->handle( name, err );
02568     MSQ_ERRZERO( err );
02569     if( !index ) MSQ_SETERR( err )( MsqError::TAG_NOT_FOUND, "could not find tag \"%s\"", name.c_str() );
02570     return (TagHandle)index;
02571 }
02572 
02573 void MeshImpl::tag_properties( TagHandle handle, std::string& name, TagType& type, unsigned& length, MsqError& err )
02574 {
02575     const TagDescription& desc = myTags->properties( (size_t)handle, err );MSQ_ERRRTN( err );
02576 
02577     name   = desc.name;
02578     type   = desc.type;
02579     length = (unsigned)( desc.size / MeshImplTags::size_from_tag_type( desc.type ) );
02580 }
02581 
02582 void MeshImpl::tag_set_element_data( TagHandle handle, size_t num_elems, const ElementHandle* elem_array,
02583                                      const void* values, MsqError& err )
02584 {
02585     myTags->set_element_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02586 }
02587 
02588 void MeshImpl::tag_get_element_data( TagHandle handle, size_t num_elems, const ElementHandle* elem_array, void* values,
02589                                      MsqError& err )
02590 {
02591     myTags->get_element_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02592 }
02593 
02594 void MeshImpl::tag_set_vertex_data( TagHandle handle, size_t num_elems, const VertexHandle* elem_array,
02595                                     const void* values, MsqError& err )
02596 {
02597     myTags->set_vertex_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02598 }
02599 
02600 void MeshImpl::tag_get_vertex_data( TagHandle handle, size_t num_elems, const VertexHandle* elem_array, void* values,
02601                                     MsqError& err )
02602 {
02603     myTags->get_vertex_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02604 }
02605 
02606 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines