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