MOAB: Mesh Oriented datABase  (version 5.3.0)
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
01177 {
01178     // deprecation of unary_function
01179     typedef size_t argument_type;
01180     typedef T result_type;
01181 
01182     T operator()( size_t idx ) const
01183     {
01184         return reinterpret_cast< T >( idx );
01185     }
01186 };
01187 
01188 void MeshImpl::vertices_get_attached_elements( const VertexHandle* vertices, size_t num_vertices,
01189                                                std::vector< ElementHandle >& elements, std::vector< size_t >& offsets,
01190                                                MsqError& err )
01191 {
01192     elements.clear();
01193     offsets.clear();
01194     size_t prev_offset = 0;
01195     offsets.reserve( num_vertices + 1 );
01196     offsets.push_back( prev_offset );
01197     const VertexHandle* const vtx_end = vertices + num_vertices;
01198     for( ; vertices < vtx_end; ++vertices )
01199     {
01200         const std::vector< size_t >& adj = myMesh->vertex_adjacencies( (size_t)*vertices, err );MSQ_ERRRTN( err );
01201 
01202         prev_offset = prev_offset + adj.size();
01203         offsets.push_back( prev_offset );
01204 
01205         std::transform( adj.begin(), adj.end(), std::back_inserter( elements ), cast_handle< ElementHandle >() );
01206     }
01207 }
01208 
01209 void MeshImpl::elements_get_attached_vertices( const ElementHandle* elements, size_t num_elems,
01210                                                std::vector< VertexHandle >& vertices, std::vector< size_t >& offsets,
01211                                                MsqError& err )
01212 {
01213     vertices.clear();
01214     offsets.clear();
01215     size_t prev_offset = 0;
01216     offsets.reserve( num_elems + 1 );
01217     offsets.push_back( prev_offset );
01218     const ElementHandle* const elem_end = elements + num_elems;
01219     for( ; elements < elem_end; ++elements )
01220     {
01221         const std::vector< size_t >& conn = myMesh->element_connectivity( (size_t)*elements, err );MSQ_ERRRTN( err );
01222 
01223         prev_offset = prev_offset + conn.size();
01224         offsets.push_back( prev_offset );
01225 
01226         std::transform( conn.begin(), conn.end(), std::back_inserter( vertices ), cast_handle< VertexHandle >() );
01227     }
01228 }
01229 
01230 // Returns the topologies of the given entities.  The "entity_topologies"
01231 // array must be at least "num_elements" in size.
01232 void MeshImpl::elements_get_topologies( const ElementHandle* element_handle_array, EntityTopology* element_topologies,
01233                                         size_t num_elements, MsqError& err )
01234 {
01235     for( size_t i = 0; i < num_elements; i++ )
01236     {
01237         element_topologies[i] = myMesh->element_topology( (size_t)element_handle_array[i], err );MSQ_CHKERR( err );
01238     }
01239 }
01240 
01241 //**************** Memory Management ****************
01242 // Tells the mesh that the client is finished with a given
01243 // entity handle.
01244 void MeshImpl::release_entity_handles( const EntityHandle* /*handle_array*/, size_t /*num_handles*/, MsqError& /*err*/ )
01245 {
01246     // Do nothing
01247 }
01248 
01249 // Instead of deleting a Mesh when you think you are done,
01250 // call release().  In simple cases, the implementation could
01251 // just call the destructor.  More sophisticated implementations
01252 // may want to keep the Mesh object to live longer than Mesquite
01253 // is using it.
01254 void MeshImpl::release()
01255 {
01256     // delete this;
01257 }
01258 
01259 const char* const vtk_type_names[] = { "bit",          "char", "unsigned_char", "short", "unsigned_short", "int",
01260                                        "unsigned_int", "long", "unsigned_long", "float", "double",         0 };
01261 
01262 void MeshImpl::read_vtk( const char* filename, MsqError& err )
01263 {
01264     int major, minor;
01265     char vendor_string[257];
01266     size_t i;
01267 
01268     FILE* file = fopen( filename, "r" );
01269     if( !file )
01270     {
01271         MSQ_SETERR( err )( MsqError::FILE_ACCESS );
01272         return;
01273     }
01274 
01275     // Read file header
01276 
01277     if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
01278     {
01279         MSQ_SETERR( err )( MsqError::IO_ERROR );
01280         fclose( file );
01281         return;
01282     }
01283 
01284     if( !strchr( vendor_string, '\n' ) || 2 != sscanf( vendor_string, "# vtk DataFile Version %d.%d", &major, &minor ) )
01285     {
01286         MSQ_SETERR( err )( MsqError::FILE_FORMAT );
01287         fclose( file );
01288         return;
01289     }
01290 
01291     if( !fgets( vendor_string, sizeof( vendor_string ), file ) )
01292     {
01293         MSQ_SETERR( err )( MsqError::IO_ERROR );
01294         fclose( file );
01295         return;
01296     }
01297 
01298     // VTK spec says this should not exceed 256 chars.
01299     if( !strchr( vendor_string, '\n' ) )
01300     {
01301         MSQ_SETERR( err )
01302         ( "Vendor string (line 2) exceeds 256 characters.", MsqError::PARSE_ERROR );
01303         fclose( file );
01304         return;
01305     }
01306 
01307     // Check file type
01308 
01309     FileTokenizer tokens( file );
01310     const char* const file_type_names[] = { "ASCII", "BINARY", 0 };
01311     int filetype                        = tokens.match_token( file_type_names, err );MSQ_ERRRTN( err );
01312     if( 2 == filetype )
01313     {
01314         MSQ_SETERR( err )
01315         ( "Cannot read BINARY VTK files -- use ASCII.", MsqError::NOT_IMPLEMENTED );
01316         return;
01317     }
01318 
01319     // Clear any existing data
01320     this->clear();
01321 
01322     const char* outer_block_names[] = { "DATASET", "FIELD", 0 };
01323     // Read the mesh
01324     // VTK docs are inconsistant with regard to whether or not
01325     // a field block should be specified as "FIELD" or "DATASET FIELD",
01326     // so allow both.
01327     while( !tokens.eof() )
01328     {
01329         int blocktype = tokens.match_token( outer_block_names, err );
01330         if( MSQ_CHKERR( err ) )
01331         {
01332             if( tokens.eof() )
01333             {
01334                 err.clear();
01335                 break;
01336             }
01337             else
01338                 return;
01339         }
01340 
01341         if( blocktype == 1 )
01342         {
01343             vtk_read_dataset( tokens, err );MSQ_ERRRTN( err );
01344         }
01345         else
01346         {
01347             vtk_read_field( tokens, err );MSQ_ERRRTN( err );
01348         }
01349     }
01350 
01351     // Make sure file actually contained some mesh
01352     if( myMesh->num_elements() == 0 )
01353     {
01354         MSQ_SETERR( err )( "File contained no mesh.", MsqError::PARSE_ERROR );
01355         return;
01356     }
01357 
01358     // There is no option for a 2-D mesh in VTK files.  Always 3
01359     numCoords = 3;
01360 
01361     // Convert tag data for fixed nodes to internal bitmap
01362     std::vector< bool > flags;
01363     tag_to_bool( "fixed", flags, err );MSQ_ERRRTN( err );
01364     if( !flags.empty() )
01365     {
01366         for( i = 0; i < myMesh->max_vertex_index(); ++i )
01367             myMesh->fix_vertex( i, flags[i], err );MSQ_ERRRTN( err );
01368     }
01369 
01370     flags.clear();
01371     tag_to_bool( "slaved", flags, err );MSQ_ERRRTN( err );
01372     if( !flags.empty() )
01373     {
01374         for( i = 0; i < myMesh->max_vertex_index(); ++i )
01375             myMesh->slave_vertex( i, flags[i], err );MSQ_ERRRTN( err );
01376     }
01377 }
01378 
01379 void MeshImpl::tag_to_bool( const char* tag_name, std::vector< bool >& values, MsqError& err )
01380 {
01381     // Convert tag data for fixed nodes to internal bitmap
01382     TagHandle handle = tag_get( tag_name, err );
01383     if( !handle || MSQ_CHKERR( err ) )
01384     {
01385         err.clear();
01386         values.clear();
01387         return;
01388     }
01389 
01390     size_t i;
01391     values.resize( myMesh->max_vertex_index(), false );
01392     const TagDescription& tag_desc = myTags->properties( (size_t)handle, err );MSQ_ERRRTN( err );
01393     bool havedata = myTags->tag_has_vertex_data( (size_t)handle, err );MSQ_ERRRTN( err );
01394     if( !havedata )
01395     {
01396         MSQ_SETERR( err )
01397         ( MsqError::FILE_FORMAT, "'%s' attribute on elements, not vertices", tag_name );
01398         return;
01399     }
01400 
01401     switch( tag_desc.type )
01402     {
01403         case BYTE: {
01404             char data;
01405             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01406             {
01407                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01408                 values[i] = !!data;
01409             }
01410             break;
01411         }
01412         case BOOL: {
01413             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01414             {
01415                 bool data;
01416                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01417                 values[i] = data;
01418             }
01419             break;
01420         }
01421         case INT: {
01422             int data;
01423             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01424             {
01425                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01426                 values[i] = !!data;
01427             }
01428             break;
01429         }
01430         case DOUBLE: {
01431             double data;
01432             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01433             {
01434                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01435                 values[i] = !!data;
01436             }
01437             break;
01438         }
01439         case HANDLE: {
01440             unsigned long data;
01441             for( i = 0; i < myMesh->max_vertex_index(); ++i )
01442             {
01443                 myTags->get_vertex_data( (size_t)handle, 1, &i, &data, err );MSQ_ERRRTN( err );
01444                 values[i] = !!data;
01445             }
01446             break;
01447         }
01448         default:
01449             MSQ_SETERR( err )( MsqError::PARSE_ERROR, "'%s' attribute has invalid type", tag_name );
01450             return;
01451     }
01452 
01453     tag_delete( handle, err );
01454 }
01455 
01456 void MeshImpl::vtk_read_dataset( FileTokenizer& tokens, MsqError& err )
01457 {
01458     const char* const data_type_names[] = {
01459         "STRUCTURED_POINTS", "STRUCTURED_GRID", "UNSTRUCTURED_GRID", "POLYDATA", "RECTILINEAR_GRID", "FIELD", 0
01460     };
01461     int datatype = tokens.match_token( data_type_names, err );MSQ_ERRRTN( err );
01462 
01463     // Ignore FIELD data at beginning of DATASET. As far as I (J.Kraftcheck)
01464     // understand the VTK documentation, there should never be a FIELD block
01465     // inside a dataset, except in attribute data.  However, some app somewhere
01466     // writes them this way, so try to handle them.
01467     for( ;; )
01468     {
01469         tokens.match_token( "FIELD", err );
01470         if( err )
01471         {
01472             tokens.unget_token();
01473             err.clear();
01474             break;
01475         }
01476         vtk_read_field( tokens, err );MSQ_ERRRTN( err );
01477     }
01478 
01479     switch( datatype )
01480     {
01481         case 1:
01482             vtk_read_structured_points( tokens, err );
01483             break;
01484         case 2:
01485             vtk_read_structured_grid( tokens, err );
01486             break;
01487         case 3:
01488             vtk_read_unstructured_grid( tokens, err );
01489             break;
01490         case 4:
01491             vtk_read_polydata( tokens, err );
01492             break;
01493         case 5:
01494             vtk_read_rectilinear_grid( tokens, err );
01495             break;
01496         case 6:
01497             vtk_read_field( tokens, err );
01498             break;
01499     }
01500 
01501     // Read attribute data
01502     const char* const block_type_names[] = { "POINT_DATA", "CELL_DATA", "DATASET", 0 };
01503     int blocktype                        = 0;
01504     while( !tokens.eof() )
01505     {
01506         // get POINT_DATA or CELL_DATA
01507         int new_block_type = tokens.match_token( block_type_names, err );
01508         if( new_block_type == 3 )  // done reading attributes
01509         {
01510             tokens.unget_token();
01511             return;
01512         }
01513 
01514         if( err )
01515         {
01516             if( tokens.eof() )
01517             {
01518                 err.clear();
01519                 break;
01520             }
01521             // If next token was neither POINT_DATA nor CELL_DATA,
01522             // then there's another attribute under the current one.
01523             if( blocktype )
01524             {
01525                 tokens.unget_token();
01526                 err.clear();
01527             }
01528             else
01529             {
01530                 MSQ_ERRRTN( err );
01531             }
01532         }
01533         else
01534         {
01535             blocktype = new_block_type;
01536             long count;
01537             tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01538 
01539             if( blocktype == 1 && (unsigned long)count != myMesh->num_vertices() )
01540             {
01541                 MSQ_SETERR( err )
01542                 ( MsqError::PARSE_ERROR,
01543                   "Count inconsistent with number of vertices"
01544                   " at line %d.",
01545                   tokens.line_number() );
01546                 return;
01547             }
01548             else if( blocktype == 2 && (unsigned long)count != myMesh->num_elements() )
01549             {
01550                 MSQ_SETERR( err )
01551                 ( MsqError::PARSE_ERROR,
01552                   "Count inconsistent with number of elements"
01553                   " at line %d.",
01554                   tokens.line_number() );
01555                 return;
01556             }
01557         }
01558 
01559         if( blocktype == 1 )
01560             vtk_read_point_data( tokens, err );
01561         else
01562             vtk_read_cell_data( tokens, err );MSQ_ERRRTN( err );
01563     }
01564 }
01565 
01566 void MeshImpl::vtk_read_structured_points( FileTokenizer& tokens, MsqError& err )
01567 {
01568     long i, j, k;
01569     long dims[3];
01570     double origin[3], space[3];
01571 
01572     tokens.match_token( "DIMENSIONS", err );MSQ_ERRRTN( err );
01573     tokens.get_long_ints( 3, dims, err );MSQ_ERRRTN( err );
01574     tokens.get_newline( err );MSQ_ERRRTN( err );
01575 
01576     if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
01577     {
01578         MSQ_SETERR( err )
01579         ( MsqError::PARSE_ERROR, "Invalid dimension at line %d", tokens.line_number() );
01580         return;
01581     }
01582 
01583     tokens.match_token( "ORIGIN", err );MSQ_ERRRTN( err );
01584     tokens.get_doubles( 3, origin, err );MSQ_ERRRTN( err );
01585     tokens.get_newline( err );MSQ_ERRRTN( err );
01586 
01587     const char* const spacing_names[] = { "SPACING", "ASPECT_RATIO", 0 };
01588     tokens.match_token( spacing_names, err );MSQ_ERRRTN( err );
01589     tokens.get_doubles( 3, space, err );MSQ_ERRRTN( err );
01590     tokens.get_newline( err );MSQ_ERRRTN( err );
01591 
01592     myMesh->allocate_vertices( dims[0] * dims[1] * dims[2], err );MSQ_ERRRTN( err );
01593     size_t vtx = 0;
01594     Vector3D off( origin[0], origin[1], origin[2] );
01595     for( k = 0; k < dims[2]; ++k )
01596         for( j = 0; j < dims[1]; ++j )
01597             for( i = 0; i < dims[0]; ++i )
01598             {
01599                 myMesh->reset_vertex( vtx++, off + Vector3D( i * space[0], j * space[1], k * space[2] ), false, err );MSQ_ERRRTN( err );
01600             }
01601 
01602     vtk_create_structured_elems( dims, err );MSQ_ERRRTN( err );
01603 }
01604 
01605 void MeshImpl::vtk_read_structured_grid( FileTokenizer& tokens, MsqError& err )
01606 {
01607     long num_verts, dims[3];
01608 
01609     tokens.match_token( "DIMENSIONS", err );MSQ_ERRRTN( err );
01610     tokens.get_long_ints( 3, dims, err );MSQ_ERRRTN( err );
01611     tokens.get_newline( err );MSQ_ERRRTN( err );
01612 
01613     if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
01614     {
01615         MSQ_SETERR( err )
01616         ( MsqError::PARSE_ERROR, "Invalid dimension at line %d", tokens.line_number() );
01617         return;
01618     }
01619 
01620     tokens.match_token( "POINTS", err );MSQ_ERRRTN( err );
01621     tokens.get_long_ints( 1, &num_verts, err );MSQ_ERRRTN( err );
01622     tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01623     tokens.get_newline( err );MSQ_ERRRTN( err );
01624 
01625     if( num_verts != ( dims[0] * dims[1] * dims[2] ) )
01626     {
01627         MSQ_SETERR( err )
01628         ( MsqError::PARSE_ERROR,
01629           "Point count not consistent with dimensions "
01630           "at line %d",
01631           tokens.line_number() );
01632         return;
01633     }
01634 
01635     myMesh->allocate_vertices( num_verts, err );MSQ_ERRRTN( err );
01636     for( size_t vtx = 0; vtx < (size_t)num_verts; ++vtx )
01637     {
01638         Vector3D pos;
01639         tokens.get_doubles( 3, const_cast< double* >( pos.to_array() ), err );MSQ_ERRRTN( err );
01640         myMesh->reset_vertex( vtx, pos, false, err );MSQ_ERRRTN( err );
01641     }
01642 
01643     vtk_create_structured_elems( dims, err );MSQ_ERRRTN( err );
01644 }
01645 
01646 void MeshImpl::vtk_read_rectilinear_grid( FileTokenizer& tokens, MsqError& err )
01647 {
01648     int i, j, k;
01649     long dims[3];
01650     const char* labels[] = { "X_COORDINATES", "Y_COORDINATES", "Z_COORDINATES" };
01651     vector< double > coords[3];
01652 
01653     tokens.match_token( "DIMENSIONS", err );MSQ_ERRRTN( err );
01654     tokens.get_long_ints( 3, dims, err );MSQ_ERRRTN( err );
01655     tokens.get_newline( err );MSQ_ERRRTN( err );
01656 
01657     if( dims[0] < 1 || dims[1] < 1 || dims[2] < 1 )
01658     {
01659         MSQ_SETERR( err )
01660         ( MsqError::PARSE_ERROR, "Invalid dimension at line %d", tokens.line_number() );
01661         return;
01662     }
01663 
01664     for( i = 0; i < 3; i++ )
01665     {
01666         long count;
01667         tokens.match_token( labels[i], err );MSQ_ERRRTN( err );
01668         tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01669         tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01670         tokens.get_newline( err );MSQ_ERRRTN( err );
01671 
01672         if( count != dims[i] )
01673         {
01674             MSQ_SETERR( err )
01675             ( MsqError::PARSE_ERROR,
01676               "Coordinate count inconsistent with dimensions"
01677               " at line %d",
01678               tokens.line_number() );
01679             return;
01680         }
01681 
01682         coords[i].resize( count );
01683         tokens.get_doubles( count, &coords[i][0], err );MSQ_ERRRTN( err );
01684     }
01685 
01686     myMesh->allocate_vertices( dims[0] * dims[1] * dims[2], err );MSQ_ERRRTN( err );
01687     size_t vtx = 0;
01688     for( k = 0; k < dims[2]; ++k )
01689         for( j = 0; j < dims[1]; ++j )
01690             for( i = 0; i < dims[0]; ++i )
01691             {
01692                 myMesh->reset_vertex( vtx++, Vector3D( coords[0][i], coords[1][j], coords[2][k] ), false, err );MSQ_ERRRTN( err );
01693             }
01694 
01695     vtk_create_structured_elems( dims, err );MSQ_ERRRTN( err );
01696 }
01697 
01698 void MeshImpl::vtk_read_polydata( FileTokenizer& tokens, MsqError& err )
01699 {
01700     long num_verts;
01701     vector< int > connectivity;
01702     const char* const poly_data_names[] = { "VERTICES", "LINES", "POLYGONS", "TRIANGLE_STRIPS", 0 };
01703 
01704     tokens.match_token( "POINTS", err );MSQ_ERRRTN( err );
01705     tokens.get_long_ints( 1, &num_verts, err );MSQ_ERRRTN( err );
01706     tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01707     tokens.get_newline( err );MSQ_ERRRTN( err );
01708 
01709     if( num_verts < 1 )
01710     {
01711         MSQ_SETERR( err )
01712         ( MsqError::PARSE_ERROR, "Invalid point count at line %d", tokens.line_number() );
01713         return;
01714     }
01715 
01716     myMesh->allocate_vertices( num_verts, err );MSQ_ERRRTN( err );
01717     for( size_t vtx = 0; vtx < (size_t)num_verts; ++vtx )
01718     {
01719         Vector3D pos;
01720         tokens.get_doubles( 3, const_cast< double* >( pos.to_array() ), err );MSQ_ERRRTN( err );
01721         myMesh->reset_vertex( vtx, pos, false, err );MSQ_ERRRTN( err );
01722     }
01723 
01724     int poly_type = tokens.match_token( poly_data_names, err );MSQ_ERRRTN( err );
01725     switch( poly_type )
01726     {
01727         case 3:
01728             vtk_read_polygons( tokens, err );MSQ_ERRRTN( err );
01729             break;
01730         case 4:
01731             MSQ_SETERR( err )
01732             ( MsqError::NOT_IMPLEMENTED, "Unsupported type: triangle strips at line %d", tokens.line_number() );
01733             return;
01734         case 1:
01735         case 2:
01736             MSQ_SETERR( err )
01737             ( MsqError::NOT_IMPLEMENTED, "Entities of dimension < 2 at line %d", tokens.line_number() );
01738             return;
01739     }
01740 }
01741 
01742 void MeshImpl::vtk_read_polygons( FileTokenizer& tokens, MsqError& err )
01743 {
01744     long size[2];
01745 
01746     tokens.get_long_ints( 2, size, err );MSQ_ERRRTN( err );
01747     tokens.get_newline( err );MSQ_ERRRTN( err );
01748     myMesh->allocate_elements( size[0], err );MSQ_ERRRTN( err );
01749     std::vector< long > conn;
01750 
01751     for( int i = 0; i < size[0]; ++i )
01752     {
01753         long count;
01754         tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01755         conn.resize( count );
01756         tokens.get_long_ints( count, arrptr( conn ), err );MSQ_ERRRTN( err );
01757         myMesh->reset_element( i, conn, POLYGON, err );MSQ_ERRRTN( err );
01758     }
01759 }
01760 
01761 void MeshImpl::vtk_read_unstructured_grid( FileTokenizer& tokens, MsqError& err )
01762 {
01763     long i, num_verts, num_elems[2];
01764 
01765     tokens.match_token( "POINTS", err );MSQ_ERRRTN( err );
01766     tokens.get_long_ints( 1, &num_verts, err );MSQ_ERRRTN( err );
01767     tokens.match_token( vtk_type_names, err );MSQ_ERRRTN( err );
01768     tokens.get_newline( err );MSQ_ERRRTN( err );
01769 
01770     if( num_verts < 1 )
01771     {
01772         MSQ_SETERR( err )
01773         ( MsqError::PARSE_ERROR, "Invalid point count at line %d", tokens.line_number() );
01774         return;
01775     }
01776 
01777     myMesh->allocate_vertices( num_verts, err );MSQ_ERRRTN( err );
01778     for( size_t vtx = 0; vtx < (size_t)num_verts; ++vtx )
01779     {
01780         Vector3D pos;
01781         tokens.get_doubles( 3, const_cast< double* >( pos.to_array() ), err );MSQ_ERRRTN( err );
01782         myMesh->reset_vertex( vtx, pos, false, err );MSQ_ERRRTN( err );
01783     }
01784 
01785     tokens.match_token( "CELLS", err );MSQ_ERRRTN( err );
01786     tokens.get_long_ints( 2, num_elems, err );MSQ_ERRRTN( err );
01787     tokens.get_newline( err );MSQ_ERRRTN( err );
01788 
01789     myMesh->allocate_elements( num_elems[0], err );MSQ_ERRRTN( err );
01790     std::vector< long > conn;
01791     for( i = 0; i < num_elems[0]; ++i )
01792     {
01793         long count;
01794         tokens.get_long_ints( 1, &count, err );MSQ_ERRRTN( err );
01795         conn.resize( count );
01796         tokens.get_long_ints( count, arrptr( conn ), err );MSQ_ERRRTN( err );
01797         myMesh->reset_element( i, conn, MIXED, err );MSQ_ERRRTN( err );
01798     }
01799 
01800     tokens.match_token( "CELL_TYPES", err );MSQ_ERRRTN( err );
01801     tokens.get_long_ints( 1, &num_elems[1], err );MSQ_ERRRTN( err );
01802     tokens.get_newline( err );MSQ_ERRRTN( err );
01803 
01804     if( num_elems[0] != num_elems[1] )
01805     {
01806         MSQ_SETERR( err )
01807         ( MsqError::PARSE_ERROR,
01808           "Number of element types does not match number of elements"
01809           "at line %d",
01810           tokens.line_number() );
01811         return;
01812     }
01813 
01814     std::vector< size_t > tconn;
01815     for( i = 0; i < num_elems[0]; ++i )
01816     {
01817         long type;
01818         size_t size;
01819         tokens.get_long_ints( 1, &type, err );MSQ_ERRRTN( err );
01820 
01821         // Check if type is a valid value
01822         const VtkTypeInfo* info = VtkTypeInfo::find_type( type, err );
01823         if( err || !info || ( !info->numNodes && type != POLYGON ) )
01824         {
01825             MSQ_SETERR( err )
01826             ( MsqError::PARSE_ERROR, "Invalid cell type %ld at line %d.", type, tokens.line_number() );
01827             return;
01828         }
01829         // Check if Mesquite supports the type
01830         if( info->msqType == MIXED )
01831         {
01832             MSQ_SETERR( err )
01833             ( MsqError::UNSUPPORTED_ELEMENT, "Unsupported cell type %ld (%s) at line %d.", type, info->name,
01834               tokens.line_number() );
01835             return;
01836         }
01837 
01838         // If node-ordering is not the same as exodus...
01839         if( info->vtkConnOrder )
01840         {
01841             size = myMesh->element_connectivity( i, err ).size();MSQ_ERRRTN( err );
01842             if( info->numNodes != size )
01843             {
01844                 MSQ_SETERR( err )
01845                 ( MsqError::UNSUPPORTED_ELEMENT, "Cell type %ld (%s) for element with %d nodes at Line %d", type,
01846                   info->name, (int)size, tokens.line_number() );
01847                 return;
01848             }
01849 
01850             tconn.resize( size );
01851             const std::vector< size_t >& pconn = myMesh->element_connectivity( i, err );MSQ_ERRRTN( err );
01852             for( size_t j = 0; j < size; ++j )
01853             {
01854                 tconn[j] = pconn[info->vtkConnOrder[j]];
01855             }
01856 
01857             myMesh->reset_element( i, tconn, info->msqType, err );MSQ_ERRRTN( err );
01858         }
01859         // Othewise (if node ordering is the same), just set the type.
01860         else
01861         {
01862             myMesh->element_topology( i, info->msqType, err );MSQ_ERRRTN( err );
01863         }
01864     }  // for(i)
01865 }
01866 
01867 void MeshImpl::vtk_create_structured_elems( const long* dims, MsqError& err )
01868 {
01869     // NOTE: this should be work fine for edges also if
01870     //      Mesquite ever supports them.  Just add the
01871     //      type for dimension 1 to the switch statement.
01872 
01873     // int non_zero[3] = {0,0,0};  // True if dim > 0 for x, y, z respectively
01874     long elem_dim  = 0;           // Element dimension (2->quad, 3->hex)
01875     long num_elems = 1;           // Total number of elements
01876     long vert_per_elem;           // Element connectivity length
01877     long edims[3] = { 1, 1, 1 };  // Number of elements in each grid direction
01878 
01879     // Populate above data
01880     for( int d = 0; d < 3; d++ )
01881         if( dims[d] > 1 )
01882         {
01883             // non_zero[elem_dim] = d;
01884             edims[d] = dims[d] - 1;
01885             num_elems *= edims[d];
01886             elem_dim++;
01887         }
01888     vert_per_elem = 1 << elem_dim;
01889 
01890     // Get element type from element dimension
01891     EntityTopology type;
01892     switch( elem_dim )
01893     {
01894             // case 1: type = EDGE;          break;
01895         case 2:
01896             type = QUADRILATERAL;
01897             break;
01898         case 3:
01899             type = HEXAHEDRON;
01900             break;
01901         default:
01902             MSQ_SETERR( err )
01903             ( "Cannot create structured mesh with elements "
01904               "of dimension < 2 or > 3.",
01905               MsqError::NOT_IMPLEMENTED );
01906             return;
01907     }
01908 
01909     // Allocate storage for elements
01910     myMesh->allocate_elements( num_elems, err );MSQ_ERRRTN( err );
01911 
01912     // Offsets of element vertices in grid relative to corner closest to origin
01913     long k                = dims[0] * dims[1];
01914     const long corners[8] = { 0, 1, 1 + dims[0], dims[0], k, k + 1, k + 1 + dims[0], k + dims[0] };
01915 
01916     // Populate element list
01917     std::vector< size_t > conn( vert_per_elem );
01918     size_t elem_idx = 0;
01919     for( long z = 0; z < edims[2]; ++z )
01920         for( long y = 0; y < edims[1]; ++y )
01921             for( long x = 0; x < edims[0]; ++x )
01922             {
01923                 const long index = x + y * dims[0] + z * ( dims[0] * dims[1] );
01924                 for( long j = 0; j < vert_per_elem; ++j )
01925                     conn[j] = index + corners[j];
01926                 myMesh->reset_element( elem_idx++, conn, type, err );MSQ_ERRRTN( err );
01927             }
01928 }
01929 
01930 void* MeshImpl::vtk_read_field_data( FileTokenizer& tokens, size_t count, size_t num_fields,
01931                                      const std::string& field_name, TagDescription& tag, MsqError& err )
01932 {
01933     tag.member = tokens.get_string( err );
01934     MSQ_ERRZERO( err );
01935     // If field is a struct containing multiple members, make
01936     // tag name the concatentation of the field name and member
01937     // name.
01938     if( num_fields > 1 )
01939     {
01940         tag.name = field_name;
01941         tag.name += " ";
01942         tag.name += tag.member;
01943         tag.member.clear();
01944     }
01945     // If field contains only one member, then make the tag
01946     // name be the field name, and store the member name to
01947     // preserve it for subsequent writes.
01948     else
01949     {
01950         tag.name = field_name;
01951     }
01952 
01953     // Get tuple size and count from the file.
01954     long sizes[2];
01955     tokens.get_long_ints( 2, sizes, err );
01956     MSQ_ERRZERO( err );
01957     if( sizes[0] < 1 )
01958     {
01959         MSQ_SETERR( err )
01960         ( MsqError::PARSE_ERROR, "Invalid tuple size (%ld) for field data %s at line %d\n", sizes[0], tag.name.c_str(),
01961           tokens.line_number() );
01962         return 0;
01963     }
01964     if( sizes[1] < 0 || ( count && (size_t)sizes[1] != count ) )
01965     {
01966         MSQ_SETERR( err )
01967         ( MsqError::PARSE_ERROR,
01968           "Invalid field data length at line %d.  "
01969           "Cannot map %lu tuples to  %ld entities.\n",
01970           tokens.line_number(), (unsigned long)count, sizes[1] );
01971         return 0;
01972     }
01973 
01974     int type = tokens.match_token( vtk_type_names, err );
01975     MSQ_ERRZERO( err );
01976     void* result = vtk_read_typed_data( tokens, type, sizes[0], sizes[1], tag, err );MSQ_CHKERR( err );
01977     return result;
01978 }
01979 
01980 void MeshImpl::vtk_read_field( FileTokenizer& tokens, MsqError& err )
01981 {
01982     // This is not supported yet.
01983     // Parse the data but throw it away because
01984     // Mesquite has no internal representation for it.
01985 
01986     std::string name = tokens.get_string( err );MSQ_ERRRTN( err );
01987     int count;
01988     tokens.get_integers( 1, &count, err );MSQ_ERRRTN( err );
01989 
01990     for( int i = 0; i < count; i++ )
01991     {
01992         TagDescription tag;
01993         void* ptr = vtk_read_field_data( tokens, 0, 1, "", tag, err );
01994         if( ptr ) free( ptr );MSQ_ERRRTN( err );
01995     }
01996 }
01997 
01998 void* MeshImpl::vtk_read_attrib_data( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
01999 {
02000     const char* const type_names[] = { "SCALARS", "COLOR_SCALARS", "VECTORS", "NORMALS", "TEXTURE_COORDINATES",
02001                                        "TENSORS", "FIELD",
02002                                        // Some buggy VTK files have empty CELL_DATA/POINT_DATA
02003                                        // blocks Try to allow them by checking for possible other
02004                                        // tokens indicating the end of the block
02005                                        "POINT_DATA", "CELL_DATA", "DATASET", 0 };
02006 
02007     int type = tokens.match_token( type_names, err );
02008     MSQ_ERRZERO( err );
02009     if( type > 7 )  // empty CELL_DATA/POINT_DATA block.
02010     {
02011         tag.vtkType = TagDescription::NONE;
02012         tokens.unget_token();
02013         return 0;
02014     }
02015 
02016     const char* name = tokens.get_string( err );
02017     MSQ_ERRZERO( err );
02018     tag.name = name;
02019 
02020     void* data = 0;
02021     switch( type )
02022     {
02023         case 1:
02024             data        = vtk_read_scalar_attrib( tokens, count, tag, err );
02025             tag.vtkType = TagDescription::SCALAR;
02026             break;
02027         case 2:
02028             data        = vtk_read_color_attrib( tokens, count, tag, err );
02029             tag.vtkType = TagDescription::COLOR;
02030             break;
02031         case 3:
02032             data        = vtk_read_vector_attrib( tokens, count, tag, err );
02033             tag.vtkType = TagDescription::VECTOR;
02034             break;
02035         case 4:
02036             data        = vtk_read_vector_attrib( tokens, count, tag, err );
02037             tag.vtkType = TagDescription::NORMAL;
02038             break;
02039         case 5:
02040             data        = vtk_read_texture_attrib( tokens, count, tag, err );
02041             tag.vtkType = TagDescription::TEXTURE;
02042             break;
02043         case 6:
02044             data        = vtk_read_tensor_attrib( tokens, count, tag, err );
02045             tag.vtkType = TagDescription::TENSOR;
02046             break;
02047         case 7:
02048             data        = 0;
02049             tag.vtkType = TagDescription::FIELD;
02050             break;
02051     }
02052 
02053     return data;
02054 }
02055 
02056 void MeshImpl::vtk_read_point_data( FileTokenizer& tokens, MsqError& err )
02057 {
02058     TagDescription tag;
02059     void* data = vtk_read_attrib_data( tokens, myMesh->num_vertices(), tag, err );
02060     if( data )  // normal attribute
02061     {
02062         vtk_store_point_data( data, tag, err );
02063         free( data );MSQ_CHKERR( err );
02064         return;
02065         ;
02066     }
02067     else if( tag.vtkType == TagDescription::FIELD )
02068     {
02069         std::string field = tag.name;
02070         int field_count;
02071         tokens.get_integers( 1, &field_count, err );MSQ_ERRRTN( err );
02072         for( int i = 0; i < field_count; ++i )
02073         {
02074             data = vtk_read_field_data( tokens, myMesh->num_vertices(), field_count, field, tag, err );MSQ_ERRRTN( err );
02075             vtk_store_point_data( data, tag, err );
02076             free( data );MSQ_ERRRTN( err );
02077         }
02078     }
02079 }
02080 
02081 void MeshImpl::vtk_store_point_data( const void* data, TagDescription& tag, MsqError& err )
02082 {
02083     size_t tag_handle = myTags->handle( tag.name, err );MSQ_ERRRTN( err );
02084     if( !tag_handle )
02085     {
02086         tag_handle = myTags->create( tag, 0, err );MSQ_ERRRTN( err );
02087     }
02088     else
02089     {
02090         const TagDescription& desc = myTags->properties( tag_handle, err );MSQ_ERRRTN( err );
02091         if( desc != tag )
02092         {
02093             MSQ_SETERR( err )
02094             ( MsqError::PARSE_ERROR,
02095               "Inconsistent types between element "
02096               "and vertex attributes named \"%s\"\n",
02097               tag.name.c_str() );
02098             return;
02099         }
02100     }
02101 
02102     std::vector< size_t > vertex_handles;
02103     myMesh->all_vertices( vertex_handles, err );MSQ_ERRRTN( err );
02104     if( !vertex_handles.empty() )
02105         myTags->set_vertex_data( tag_handle, vertex_handles.size(), arrptr( vertex_handles ), data, err );MSQ_ERRRTN( err );
02106 }
02107 
02108 void MeshImpl::vtk_read_cell_data( FileTokenizer& tokens, MsqError& err )
02109 {
02110     TagDescription tag;
02111     void* data = vtk_read_attrib_data( tokens, myMesh->num_elements(), tag, err );
02112     if( data )  // normal attribute
02113     {
02114         vtk_store_cell_data( data, tag, err );
02115         free( data );MSQ_CHKERR( err );
02116         return;
02117         ;
02118     }
02119     else if( tag.vtkType == TagDescription::FIELD )
02120     {
02121         std::string field = tag.name;
02122         int field_count;
02123         tokens.get_integers( 1, &field_count, err );MSQ_ERRRTN( err );
02124         for( int i = 0; i < field_count; ++i )
02125         {
02126             data = vtk_read_field_data( tokens, myMesh->num_elements(), field_count, field, tag, err );MSQ_ERRRTN( err );
02127             vtk_store_cell_data( data, tag, err );
02128             free( data );MSQ_ERRRTN( err );
02129         }
02130     }
02131 }
02132 
02133 void MeshImpl::vtk_store_cell_data( const void* data, TagDescription& tag, MsqError& err )
02134 {
02135     size_t tag_handle = myTags->handle( tag.name, err );MSQ_ERRRTN( err );
02136     if( !tag_handle )
02137     {
02138         tag_handle = myTags->create( tag, 0, err );MSQ_ERRRTN( err );
02139     }
02140     else
02141     {
02142         const TagDescription& desc = myTags->properties( tag_handle, err );MSQ_ERRRTN( err );
02143         if( desc != tag )
02144         {
02145             MSQ_SETERR( err )
02146             ( MsqError::PARSE_ERROR,
02147               "Inconsistent types between element "
02148               "and vertex attributes named \"%s\"\n",
02149               tag.name.c_str() );
02150             return;
02151         }
02152     }
02153 
02154     std::vector< size_t > element_handles;
02155     myMesh->all_elements( element_handles, err );MSQ_ERRRTN( err );
02156     if( !element_handles.empty() )
02157         myTags->set_element_data( tag_handle, element_handles.size(), arrptr( element_handles ), data, err );MSQ_ERRRTN( err );
02158 }
02159 
02160 void* MeshImpl::vtk_read_typed_data( FileTokenizer& tokens, int type, size_t per_elem, size_t num_elem,
02161                                      TagDescription& tag, MsqError& err )
02162 {
02163     void* data_ptr;
02164     size_t count = per_elem * num_elem;
02165     switch( type )
02166     {
02167         case 1:
02168             tag.size = per_elem * sizeof( bool );
02169             tag.type = BOOL;
02170             data_ptr = malloc( num_elem * tag.size );
02171             tokens.get_booleans( count, (bool*)data_ptr, err );
02172             break;
02173         case 2:
02174         case 3:
02175         case 4:
02176         case 5:
02177         case 6:
02178         case 7:
02179             tag.size = per_elem * sizeof( int );
02180             tag.type = INT;
02181             data_ptr = malloc( num_elem * tag.size );
02182             tokens.get_integers( count, (int*)data_ptr, err );
02183             break;
02184         case 8:
02185         case 9:
02186             // this is a bit of a hack since MeshImpl doesn't have a LONG type (HANDLE is used by
02187             // ParallelMesh for long)
02188             tag.size = per_elem * sizeof( size_t );
02189             assert( sizeof( long ) == sizeof( size_t ) );
02190             assert( sizeof( long ) == sizeof( void* ) );
02191             tag.type = HANDLE;
02192             data_ptr = malloc( num_elem * tag.size );
02193             tokens.get_long_ints( count, (long*)data_ptr, err );
02194             break;
02195         case 10:
02196         case 11:
02197             tag.size = per_elem * sizeof( double );
02198             tag.type = DOUBLE;
02199             data_ptr = malloc( num_elem * tag.size );
02200             tokens.get_doubles( count, (double*)data_ptr, err );
02201             break;
02202         default:
02203             MSQ_SETERR( err )( "Invalid data type", MsqError::INVALID_ARG );
02204             return 0;
02205     }
02206 
02207     if( MSQ_CHKERR( err ) )
02208     {
02209         free( data_ptr );
02210         return 0;
02211     }
02212 
02213     return data_ptr;
02214 }
02215 
02216 void* MeshImpl::vtk_read_scalar_attrib( FileTokenizer& tokens, long count, TagDescription& desc, MsqError& err )
02217 {
02218     if( !count ) return 0;
02219 
02220     int type = tokens.match_token( vtk_type_names, err );
02221     MSQ_ERRZERO( err );
02222 
02223     long size;
02224     const char* tok = tokens.get_string( err );
02225     MSQ_ERRZERO( err );
02226     const char* end = 0;
02227     size            = strtol( tok, (char**)&end, 0 );
02228     if( *end )
02229     {
02230         size = 1;
02231         tokens.unget_token();
02232     }
02233 
02234     // VTK spec says cannot be greater than 4--do we care?
02235     if( size < 1 || size > 4 )
02236     {
02237         MSQ_SETERR( err )
02238         ( MsqError::PARSE_ERROR,
02239           "Scalar count out of range [1,4]"
02240           " at line %d",
02241           tokens.line_number() );
02242         return 0;
02243     }
02244 
02245     tokens.match_token( "LOOKUP_TABLE", err );
02246     MSQ_ERRZERO( err );
02247     tok = tokens.get_string( err );
02248     MSQ_ERRZERO( err );
02249 
02250     // If no lookup table, just read and return the data
02251     if( !strcmp( tok, "default" ) )
02252     {
02253         void* ptr = vtk_read_typed_data( tokens, type, size, count, desc, err );
02254         MSQ_ERRZERO( err );
02255         return ptr;
02256     }
02257 
02258     // If we got this far, then the data has a lookup
02259     // table.  First read the lookup table and convert
02260     // to integers.
02261     string name = tok;
02262     vector< long > table( size * count );
02263     if( type > 0 && type < 10 )  // Is an integer type
02264     {
02265         tokens.get_long_ints( table.size(), arrptr( table ), err );
02266         MSQ_ERRZERO( err );
02267     }
02268     else  // Is a real-number type
02269     {
02270         for( std::vector< long >::iterator iter = table.begin(); iter != table.end(); ++iter )
02271         {
02272             double data;
02273             tokens.get_doubles( 1, &data, err );
02274             MSQ_ERRZERO( err );
02275 
02276             *iter = (long)data;
02277             if( (double)*iter != data )
02278             {
02279                 MSQ_SETERR( err )
02280                 ( MsqError::PARSE_ERROR, "Invalid lookup index (%.0f) at line %d", data, tokens.line_number() );
02281                 return 0;
02282             }
02283         }
02284     }
02285 
02286     // Now read the data - must be float RGBA color triples
02287 
02288     long table_count;
02289     tokens.match_token( "LOOKUP_TABLE", err );
02290     MSQ_ERRZERO( err );
02291     tokens.match_token( name.c_str(), err );
02292     MSQ_ERRZERO( err );
02293     tokens.get_long_ints( 1, &table_count, err );
02294     MSQ_ERRZERO( err );
02295 
02296     vector< float > table_data( table_count * 4 );
02297     tokens.get_floats( table_data.size(), arrptr( table_data ), err );
02298     MSQ_ERRZERO( err );
02299 
02300     // Create list from indexed data
02301 
02302     float* data      = (float*)malloc( sizeof( float ) * count * size * 4 );
02303     float* data_iter = data;
02304     for( std::vector< long >::iterator idx = table.begin(); idx != table.end(); ++idx )
02305     {
02306         if( *idx < 0 || *idx >= table_count )
02307         {
02308             MSQ_SETERR( err )
02309             ( MsqError::PARSE_ERROR, "LOOKUP_TABLE index %ld out of range.", *idx );
02310             free( data );
02311             return 0;
02312         }
02313 
02314         for( int i = 0; i < 4; i++ )
02315         {
02316             *data_iter = table_data[4 * *idx + i];
02317             ++data_iter;
02318         }
02319     }
02320 
02321     desc.size = size * 4 * sizeof( float );
02322     desc.type = DOUBLE;
02323     return data;
02324 }
02325 
02326 void* MeshImpl::vtk_read_color_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02327 {
02328     long size;
02329     tokens.get_long_ints( 1, &size, err );
02330     MSQ_ERRZERO( err );
02331 
02332     if( size < 1 )
02333     {
02334         MSQ_SETERR( err )
02335         ( MsqError::PARSE_ERROR, "Invalid size (%ld) at line %d", size, tokens.line_number() );
02336         return 0;
02337     }
02338 
02339     float* data = (float*)malloc( sizeof( float ) * count * size );
02340     tokens.get_floats( count * size, data, err );
02341     if( MSQ_CHKERR( err ) )
02342     {
02343         free( data );
02344         return 0;
02345     }
02346 
02347     tag.size = size * sizeof( float );
02348     tag.type = DOUBLE;
02349     return data;
02350 }
02351 
02352 void* MeshImpl::vtk_read_vector_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02353 {
02354     int type = tokens.match_token( vtk_type_names, err );
02355     MSQ_ERRZERO( err );
02356 
02357     void* result = vtk_read_typed_data( tokens, type, 3, count, tag, err );
02358     MSQ_ERRZERO( err );
02359     return result;
02360 }
02361 
02362 void* MeshImpl::vtk_read_texture_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02363 {
02364     int type, dim;
02365     tokens.get_integers( 1, &dim, err );
02366     MSQ_ERRZERO( err );
02367     type = tokens.match_token( vtk_type_names, err );
02368     MSQ_ERRZERO( err );
02369 
02370     if( dim < 1 || dim > 3 )
02371     {
02372         MSQ_SETERR( err )
02373         ( MsqError::PARSE_ERROR, "Invalid dimension (%d) at line %d.", dim, tokens.line_number() );
02374         return 0;
02375     }
02376 
02377     void* result = vtk_read_typed_data( tokens, type, dim, count, tag, err );
02378     MSQ_ERRZERO( err );
02379     return result;
02380 }
02381 
02382 void* MeshImpl::vtk_read_tensor_attrib( FileTokenizer& tokens, long count, TagDescription& tag, MsqError& err )
02383 {
02384     int type = tokens.match_token( vtk_type_names, err );
02385     MSQ_ERRZERO( err );
02386 
02387     void* result = vtk_read_typed_data( tokens, type, 9, count, tag, err );
02388     MSQ_ERRZERO( err );
02389     return result;
02390 }
02391 
02392 void MeshImpl::vtk_write_attrib_data( std::ostream& file, const TagDescription& desc, const void* data, size_t count,
02393                                       MsqError& err ) const
02394 {
02395     // srkenno@sandia.gov: we now allow this type to be able to write e.g. GLOBAL_ID for parallel
02396     // meshes
02397     /*
02398     if (desc.type == HANDLE)
02399     {
02400       MSQ_SETERR(err)("Cannot write HANDLE tag data to VTK file.",
02401                       MsqError::FILE_FORMAT);
02402       return;
02403       }*/
02404 
02405     TagDescription::VtkType vtk_type = desc.vtkType;
02406     unsigned vlen                    = desc.size / MeshImplTags::size_from_tag_type( desc.type );
02407     // guess one from data length if not set
02408     if( vtk_type == TagDescription::NONE )
02409     {
02410         switch( vlen )
02411         {
02412             default:
02413                 vtk_type = TagDescription::SCALAR;
02414                 break;
02415             case 3:
02416                 vtk_type = TagDescription::VECTOR;
02417                 break;
02418             case 9:
02419                 vtk_type = TagDescription::TENSOR;
02420                 break;
02421                 return;
02422         }
02423     }
02424 
02425     // srkenno@sandia.gov: from class Mesh, the typenames below should correspond in order...
02426     // enum TagType { BYTE, BOOL, INT, DOUBLE, HANDLE };
02427 
02428     const char* const typenames[] = { "unsigned_char", "bit", "int", "double", "unsigned_long" };
02429     std::string field, member;
02430 
02431     int num_per_line;
02432     switch( vtk_type )
02433     {
02434         case TagDescription::SCALAR:
02435             num_per_line = vlen;
02436             file << "SCALARS " << desc.name << " " << typenames[desc.type] << " " << vlen << "\n";
02437             file << "LOOKUP_TABLE default\n";
02438             break;
02439         case TagDescription::COLOR:
02440             num_per_line = vlen;
02441             file << "COLOR_SCALARS " << desc.name << " " << vlen << "\n";
02442             break;
02443         case TagDescription::VECTOR:
02444             num_per_line = 3;
02445             if( vlen != 3 )
02446             {
02447                 MSQ_SETERR( err )
02448                 ( MsqError::INTERNAL_ERROR, "Tag \"%s\" is labeled as a VTK vector attribute but has %u values.",
02449                   desc.name.c_str(), vlen );
02450                 return;
02451             }
02452             file << "VECTORS " << desc.name << " " << typenames[desc.type] << "\n";
02453             break;
02454         case TagDescription::NORMAL:
02455             num_per_line = 3;
02456             if( vlen != 3 )
02457             {
02458                 MSQ_SETERR( err )
02459                 ( MsqError::INTERNAL_ERROR, "Tag \"%s\" is labeled as a VTK normal attribute but has %u values.",
02460                   desc.name.c_str(), vlen );
02461                 return;
02462             }
02463             file << "NORMALS " << desc.name << " " << typenames[desc.type] << "\n";
02464             break;
02465         case TagDescription::TEXTURE:
02466             num_per_line = vlen;
02467             file << "TEXTURE_COORDINATES " << desc.name << " " << typenames[desc.type] << " " << vlen << "\n";
02468             break;
02469         case TagDescription::TENSOR:
02470             num_per_line = 3;
02471             if( vlen != 9 )
02472             {
02473                 MSQ_SETERR( err )
02474                 ( MsqError::INTERNAL_ERROR, "Tag \"%s\" is labeled as a VTK tensor attribute but has %u values.",
02475                   desc.name.c_str(), vlen );
02476                 return;
02477             }
02478             file << "TENSORS " << desc.name << " " << typenames[desc.type] << "\n";
02479             break;
02480         case TagDescription::FIELD:
02481             num_per_line = vlen;
02482             get_field_names( desc, field, member, err );MSQ_ERRRTN( err );
02483             file << member << " " << vlen << " " << count << " " << typenames[desc.type] << "\n";
02484             break;
02485         default:
02486             MSQ_SETERR( err )( "Unknown VTK attribute type for tag.", MsqError::INTERNAL_ERROR );
02487             return;
02488     }
02489 
02490     size_t i, total = count * vlen;
02491     char* space = new char[num_per_line];
02492     memset( space, ' ', num_per_line );
02493     space[num_per_line - 1]    = '\n';
02494     const unsigned char* odata = (const unsigned char*)data;
02495     const bool* bdata          = (const bool*)data;
02496     const int* idata           = (const int*)data;
02497     const long* ldata          = (const long*)data;
02498     const double* ddata        = (const double*)data;
02499     switch( desc.type )
02500     {
02501         case BYTE:
02502             for( i = 0; i < total; ++i )
02503                 file << (unsigned int)odata[i] << space[i % num_per_line];
02504             break;
02505         case BOOL:
02506             for( i = 0; i < total; ++i )
02507                 file << ( bdata[i] ? '1' : '0' ) << space[i % num_per_line];
02508             break;
02509         case INT:
02510             for( i = 0; i < total; ++i )
02511                 file << idata[i] << space[i % num_per_line];
02512             break;
02513         case DOUBLE:
02514             for( i = 0; i < total; ++i )
02515                 file << ddata[i] << space[i % num_per_line];
02516             break;
02517         case HANDLE:
02518             for( i = 0; i < total; ++i )
02519                 file << ldata[i] << space[i % num_per_line];
02520             break;
02521         default:
02522             MSQ_SETERR( err )( "Unknown tag type.", MsqError::INTERNAL_ERROR );
02523     }
02524     delete[] space;
02525 }
02526 
02527 /**************************************************************************
02528  *                               TAGS
02529  **************************************************************************/
02530 
02531 TagHandle MeshImpl::tag_create( const std::string& name, TagType type, unsigned length, const void* defval,
02532                                 MsqError& err )
02533 {
02534     TagDescription::VtkType vtype;
02535     std::string field;
02536     switch( length )
02537     {
02538         case 1:
02539             vtype = TagDescription::SCALAR;
02540             break;
02541         case 3:
02542             vtype = TagDescription::VECTOR;
02543             break;
02544         case 9:
02545             vtype = TagDescription::TENSOR;
02546             break;
02547         default:
02548             vtype = TagDescription::FIELD;
02549             field = MESQUITE_FIELD_TAG;
02550             break;
02551     }
02552 
02553     // If tag name contains a space, assume the tag name
02554     // is a concatenation of the VTK field and member names.
02555     if( vtype != TagDescription::FIELD && name.find( " " ) != std::string::npos ) vtype = TagDescription::FIELD;
02556 
02557     size_t size = MeshImplTags::size_from_tag_type( type );
02558     TagDescription desc( name, type, vtype, length * size, field );
02559     size_t index = myTags->create( desc, defval, err );
02560     MSQ_ERRZERO( err );
02561     return (TagHandle)index;
02562 }
02563 
02564 void MeshImpl::tag_delete( TagHandle handle, MsqError& err )
02565 {
02566     myTags->destroy( (size_t)handle, err );MSQ_CHKERR( err );
02567 }
02568 
02569 TagHandle MeshImpl::tag_get( const std::string& name, MsqError& err )
02570 {
02571     size_t index = myTags->handle( name, err );
02572     MSQ_ERRZERO( err );
02573     if( !index ) MSQ_SETERR( err )( MsqError::TAG_NOT_FOUND, "could not find tag \"%s\"", name.c_str() );
02574     return (TagHandle)index;
02575 }
02576 
02577 void MeshImpl::tag_properties( TagHandle handle, std::string& name, TagType& type, unsigned& length, MsqError& err )
02578 {
02579     const TagDescription& desc = myTags->properties( (size_t)handle, err );MSQ_ERRRTN( err );
02580 
02581     name   = desc.name;
02582     type   = desc.type;
02583     length = (unsigned)( desc.size / MeshImplTags::size_from_tag_type( desc.type ) );
02584 }
02585 
02586 void MeshImpl::tag_set_element_data( TagHandle handle, size_t num_elems, const ElementHandle* elem_array,
02587                                      const void* values, MsqError& err )
02588 {
02589     myTags->set_element_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02590 }
02591 
02592 void MeshImpl::tag_get_element_data( TagHandle handle, size_t num_elems, const ElementHandle* elem_array, void* values,
02593                                      MsqError& err )
02594 {
02595     myTags->get_element_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02596 }
02597 
02598 void MeshImpl::tag_set_vertex_data( TagHandle handle, size_t num_elems, const VertexHandle* elem_array,
02599                                     const void* values, MsqError& err )
02600 {
02601     myTags->set_vertex_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02602 }
02603 
02604 void MeshImpl::tag_get_vertex_data( TagHandle handle, size_t num_elems, const VertexHandle* elem_array, void* values,
02605                                     MsqError& err )
02606 {
02607     myTags->get_vertex_data( (size_t)handle, num_elems, (const size_t*)elem_array, values, err );MSQ_CHKERR( err );
02608 }
02609 
02610 }  // namespace MBMesquite
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines