C++ MeshImpl::read_vtk方法代码示例

    // Test reading Vtk vector attribute
  void VtkTest::test_read_vector_attrib()
  {
    MeshImpl mesh;
    MsqPrintError err(cout);
    
    FILE* file = fopen( temp_file_name, "w+" );
    fputs( structured_3d_points_data, file );
    fputs( simple_vector_attrib, file );
    fclose( file );
    
    mesh.read_vtk( temp_file_name, err );
    remove( temp_file_name );
    ASSERT_NO_ERROR(err);
    
    std::vector<Mesh::ElementHandle> elems;
    mesh.get_all_elements( elems, err );
    CPPUNIT_ASSERT( !err );
    CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
   
    void* th = mesh.tag_get( "hexvect", err );
    CPPUNIT_ASSERT( !err );

    std::string name;
    Mesh::TagType type;
    unsigned tagsize;
    mesh.tag_properties( th, name, type, tagsize, err );
    CPPUNIT_ASSERT( !err && type == Mesh::DOUBLE && tagsize == 3 );
    
    double elem_data[24];
    mesh.tag_get_element_data( th, 8, arrptr(elems), elem_data, err );
    CPPUNIT_ASSERT( !err );
    
    for (int i = 0; i < 8; ++i)
      CPPUNIT_ASSERT( Vector3D( elem_data+3*i ) == Vector3D( i+1, i+1, i+1 ) );
  }

  void VtkTest::check_field_attrib( const char* temp_file_name )
  {
    MeshImpl mesh;
    MsqPrintError err(cout);
    
    mesh.read_vtk( temp_file_name, err );
    remove( temp_file_name );
    ASSERT_NO_ERROR(err);
    
    std::vector<Mesh::ElementHandle> elems;
    mesh.get_all_elements( elems, err );
    CPPUNIT_ASSERT( !err );
    CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
   
    std::string name;
    Mesh::TagType type;
    unsigned tagsize;


    void* th = mesh.tag_get( "test_field elem_vects", err );
    CPPUNIT_ASSERT( !err );

    mesh.tag_properties( th, name, type, tagsize, err );
    CPPUNIT_ASSERT( !err && type == Mesh::DOUBLE && tagsize == 3 );
    
    double elem_data[24];
    mesh.tag_get_element_data( th, 8, arrptr(elems), elem_data, err );
    CPPUNIT_ASSERT( !err );
    
    for (int i = 0; i < 8; ++i)
      CPPUNIT_ASSERT( Vector3D( elem_data+3*i ) == Vector3D( i+1, i+1, i+1 ) );
   
    
    th = mesh.tag_get( "test_field elem_ids", err );
    CPPUNIT_ASSERT( !err );

    mesh.tag_properties( th, name, type, tagsize, err );
    CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
    
    int elem_ids[8];
    mesh.tag_get_element_data( th, 8, arrptr(elems), elem_ids, err );
    CPPUNIT_ASSERT( !err );
    
    for (int i = 0; i < 8; ++i)
      CPPUNIT_ASSERT( elem_ids[i] == i+1 );
   
    
    th = mesh.tag_get( "field1", err );
    CPPUNIT_ASSERT( !err );

    mesh.tag_properties( th, name, type, tagsize, err );
    CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
    
    int values[8];
    mesh.tag_get_element_data( th, 8, arrptr(elems), values, err );
    CPPUNIT_ASSERT( !err );
    
    for (int i = 0; i < 8; ++i)
      CPPUNIT_ASSERT( values[i] == 8-i );
  }

void TagVertexMeshTest::setUp()
{
  const char vtk_data[] = 
    "# vtk DataFile Version 2.0\n"
    "test mesh\n"
    "ASCII\n"
    "DATASET UNSTRUCTURED_GRID\n"
    "POINTS 3 float\n"
    "0 0 0\n"
    "1 0 0\n"
    "0 1 0\n"
    "CELLS 1 4\n"
    "3 0 1 2\n"
    "CELL_TYPES 1\n"
    "5\n";
  
  FILE* file = fopen( TEMP_FILE_NAME, "w" );
  CPPUNIT_ASSERT( !!file );
  size_t r = fwrite( vtk_data, sizeof(vtk_data)-1, 1, file );
  fclose( file );
  CPPUNIT_ASSERT( r == 1 );
  
  MsqPrintError err( std::cerr );
  realMesh = new MeshImpl;
  realMesh->read_vtk( TEMP_FILE_NAME, err );
  remove( TEMP_FILE_NAME );
  ASSERT_NO_ERROR(err);
}

void SlaveBoundaryVerticesTest::make_mesh( MeshImpl& mesh,
                                           DomainClassifier& domain,
                                           const int intervals )
{
  MsqPrintError err(std::cerr);
  const char input_file[] = MESH_FILES_DIR "3D/VTK/6x6x6-hex20.vtk";
  
  const Vector3D min( -3, -3, -3 );
  const Vector3D max(  3,  3,  3 );
  const Vector3D space = (max - min) * 1.0/(intervals);
  
  mesh.clear();
  mesh.read_vtk( input_file, err );
  ASSERT_NO_ERROR(err);
  
  const Vector3D corners[8] = { Vector3D(min[0],min[1],min[2]),
                                Vector3D(max[0],min[1],min[2]),
                                Vector3D(max[0],max[1],min[2]),
                                Vector3D(min[0],max[1],min[2]),
                                Vector3D(min[0],min[1],max[2]),
                                Vector3D(max[0],min[1],max[2]),
                                Vector3D(max[0],max[1],max[2]),
                                Vector3D(min[0],max[1],max[2]) };
  
  MeshDomain* subdomains[26] = { 
                    new PlanarDomain( PlanarDomain::XZ, min[1] ),
                    new PlanarDomain( PlanarDomain::YZ, max[0] ),
                    new PlanarDomain( PlanarDomain::XZ, max[1] ),
                    new PlanarDomain( PlanarDomain::YZ, min[0] ),
                    new PlanarDomain( PlanarDomain::XY, min[2] ),
                    new PlanarDomain( PlanarDomain::XY, max[2] ),
                    new LineDomain( corners[0], corners[1] - corners[0] ),
                    new LineDomain( corners[1], corners[2] - corners[1] ),
                    new LineDomain( corners[2], corners[3] - corners[2] ),
                    new LineDomain( corners[3], corners[0] - corners[3] ),
                    new LineDomain( corners[0], corners[4] - corners[0] ),
                    new LineDomain( corners[1], corners[5] - corners[1] ),
                    new LineDomain( corners[2], corners[6] - corners[0] ),
                    new LineDomain( corners[3], corners[7] - corners[1] ),
                    new LineDomain( corners[4], corners[5] - corners[4] ),
                    new LineDomain( corners[5], corners[6] - corners[5] ),
                    new LineDomain( corners[6], corners[7] - corners[6] ),
                    new LineDomain( corners[7], corners[4] - corners[7] ),
                    new PointDomain( corners[0] ),
                    new PointDomain( corners[1] ),
                    new PointDomain( corners[2] ),
                    new PointDomain( corners[3] ),
                    new PointDomain( corners[4] ),
                    new PointDomain( corners[5] ),
                    new PointDomain( corners[6] ),
                    new PointDomain( corners[7] ) };
  const int subdims[26] = { 2, 2, 2, 2, 2, 2,
                            1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                            0, 0, 0, 0, 0, 0, 0, 0 };
  DomainClassifier::classify_skin_geometrically( domain, &mesh, 1e-6,
                                                 subdomains, subdims, 26,
                                                 err );
  domain.delete_sub_domains(true);
  ASSERT_NO_ERROR(err);
}

    // Test reading Vtk simple (one-component) scalar attribute
  void VtkTest::test_read_simple_scalar_attrib()
  {
    MeshImpl mesh;
    MsqPrintError err(cout);
   
    FILE* file = fopen( temp_file_name, "w+" );
    fputs( structured_3d_points_data, file );
    fputs( simple_scalar_attrib, file );
    fclose( file );
    
    mesh.read_vtk( temp_file_name, err );
    remove( temp_file_name );
    ASSERT_NO_ERROR(err);
    
    std::vector<Mesh::ElementHandle> elems;
    mesh.get_all_elements( elems, err );
    CPPUNIT_ASSERT( !err );
    CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
   
    void* th = mesh.tag_get( "global_id", err );
    CPPUNIT_ASSERT( !err );

    std::string name;
    Mesh::TagType type;
    unsigned tagsize;
    mesh.tag_properties( th, name, type, tagsize, err );
    CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
    
    int elem_data[8];
    mesh.tag_get_element_data( th, 8, arrptr(elems), elem_data, err );
    CPPUNIT_ASSERT( !err );
    
    for (int i = 0; i < 8; ++i)
      CPPUNIT_ASSERT( elem_data[i] == (1+i) );
  }

   // Test writing VTK unstructured mesh
 void VtkTest::test_write()
 {
   MeshImpl mesh;
   MsqPrintError err(cout);
   
     // Create file containing unstructured mesh test case
   FILE* file = fopen( temp_file_name, "w+" );
   CPPUNIT_ASSERT(file);
   int rval = fputs( mixed_unstructured_data, file );
   fclose( file );
   if (rval == EOF) remove(temp_file_name);
   CPPUNIT_ASSERT(rval != EOF);
   
     // Read unstructured mesh file
   mesh.read_vtk( temp_file_name, err );
   remove(temp_file_name);
   ASSERT_NO_ERROR(err);
   
     // Write unstructured mesh file back out
   mesh.write_vtk( temp_file_name, err );
   if (err)  remove( temp_file_name );
   ASSERT_NO_ERROR(err);
   
     // Check if file contained expected mesh
   test_read_unstructured( temp_file_name );
   remove( temp_file_name );
 }

void get_planar_example( const char* filename,
                         DomainClassifier& geom,
                         MeshImpl& mesh,
                         MsqError& err )
{
  static PlanarDomain z(PlanarDomain::XY);
  mesh.read_vtk(filename, err); MSQ_ERRRTN(err);
  mesh.mark_skin_fixed(err); MSQ_ERRRTN(err);
  DomainClassifier::DomainSet set(&z);
  mesh.get_all_vertices( set.vertices, err ); MSQ_ERRRTN(err);
  mesh.get_all_elements( set.elements, err ); MSQ_ERRRTN(err);
  DomainClassifier::classify_by_handle( geom, &mesh, &set, 1, err ); MSQ_ERRRTN(err);
}

   // Test reading 3D Vtk rectilinear-grid mesh
 void VtkTest::test_read_rectilinear_grid()
 {
   MeshImpl mesh;
   MsqPrintError err(cout);
   
   FILE* file = fopen( temp_file_name, "w+" );
   fputs( rectilinear_grid_data, file );
   fclose( file );
   
   mesh.read_vtk( temp_file_name, err );
   remove( temp_file_name );
   ASSERT_NO_ERROR(err);
   
   check_8hex_structured( mesh );
 }

   // Test reading 2D Vtk structured-points mesh
 void VtkTest::test_read_structured_2d_points()
 {
   MeshImpl mesh;
   MsqPrintError err(cout);
  
   FILE* file = fopen( temp_file_name, "w+" );
   fputs( structured_2d_points_data, file );
   fclose( file );
   
   mesh.read_vtk( temp_file_name, err );
   remove( temp_file_name );
   ASSERT_NO_ERROR(err);
   
   check_4quad_structured( mesh );
 }

  void VtkTest::test_elements()
  {
    Mesquite::MsqPrintError err(cout);
    MeshImpl mMesh;
    mMesh.read_vtk(MESH_FILES_DIR "2D/vtk/tris/untangled/equil_tri2.vtk", err);
    ASSERT_NO_ERROR(err);
    Mesquite::MeshDomainAssoc mesh_and_domain = Mesquite::MeshDomainAssoc(&mMesh, 0);
    Mesquite::Instruction::initialize_vertex_byte( &mesh_and_domain, 0, err );
    ASSERT_NO_ERROR(err);
    
      // Retrieve a patch
    Mesquite::PatchData pd;
    pd.set_mesh( &mMesh );
    VertexPatches patch_set;
    patch_set.set_mesh( &mMesh );
    PatchIterator patches( &patch_set );
    patches.get_next_patch( pd, err );
    ASSERT_NO_ERROR(err);
    
    int free_vtx = pd.num_free_vertices(); 
//    std::cout << "nb of free vertices: " << free_vtx << std::endl;
    CPPUNIT_ASSERT( free_vtx == 1 );
    
    Mesquite::MsqMeshEntity* element_array =  pd.get_element_array(err); ASSERT_NO_ERROR(err);
    size_t num_elements = pd.num_elements();
    CPPUNIT_ASSERT( num_elements == 6 );
    
    const Mesquite::MsqVertex* vtx_array = pd.get_vertex_array(err); ASSERT_NO_ERROR(err);
    size_t num_vertices = pd.num_nodes();
    CPPUNIT_ASSERT( num_vertices == 7 );
    
    CPPUNIT_ASSERT( tri_check_validity(element_array, num_elements, vtx_array, num_vertices) == 1 );
    
    patches.get_next_patch( pd, err ); ASSERT_NO_ERROR(err);
    
    element_array =  pd.get_element_array(err); ASSERT_NO_ERROR(err);
    num_elements = pd.num_elements();
    CPPUNIT_ASSERT( num_elements == 6 );
    
    vtx_array = pd.get_vertex_array(err); ASSERT_NO_ERROR(err);
    num_vertices = pd.num_nodes();
    CPPUNIT_ASSERT( num_vertices == 7 );
    
    CPPUNIT_ASSERT( tri_check_validity(element_array, num_elements, vtx_array, num_vertices) == 1 );
  }

    // Test reading MeshImpl boundary-vertex bit
    // from Vtk scalar attribute.
  void VtkTest::test_read_fixed_attrib()
  {
    MeshImpl mesh;
    MsqPrintError err(cout);
    
    FILE* file = fopen( temp_file_name, "w+" );
    fputs( structured_3d_points_data, file );
    fputs( fixed_vertex_attrib, file );
    fclose( file );
    
    mesh.read_vtk( temp_file_name, err );
    remove( temp_file_name );
    ASSERT_NO_ERROR(err);

    std::vector<Mesh::ElementHandle> elems;
    mesh.get_all_elements( elems, err );
    CPPUNIT_ASSERT( !err );
    CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
    
    std::vector<Mesh::VertexHandle> verts;
    std::vector<size_t> offsets;
    mesh.elements_get_attached_vertices( arrptr(elems), elems.size(), verts, offsets, err );
    ASSERT_NO_ERROR(err);
    
    // get unique list of vertices
    std::vector<Mesh::VertexHandle>::iterator new_end;
    std::sort( verts.begin(), verts.end() );
    new_end = std::unique( verts.begin(), verts.end() );
    verts.resize( new_end - verts.begin() );
    CPPUNIT_ASSERT_EQUAL( verts.size(), (size_t)27 );

    // get fixed flag
    std::vector<bool> fixed;
    mesh.vertices_get_fixed_flag( arrptr(verts), fixed, verts.size(), err );
    ASSERT_NO_ERROR(err);
    CPPUNIT_ASSERT_EQUAL( verts.size(), fixed.size() );
    
    for (int i = 0; i < 27; ++i)
    {
      bool should_be_fixed = (i != 13);
      CPPUNIT_ASSERT_EQUAL( should_be_fixed, (bool)fixed[i] );
    }
  }

std::string get_sphere_cylinder_example( DomainClassifier& geom, MeshImpl& mesh, MsqError& err )
{
  const char filename[] = SRCDIR "sphereCylinder_1194_inv.vtk";

  const Vector3D vec_k(0,0,8), vec_nk(0,0,-8);
  const Vector3D vec_c(0,0,5);

  //++ 0D domains ++

  //++ 1D domains ++

  //top circle
  static CircleDomain cr_to(vec_k,vec_k,8.0), cr_ti(vec_k,vec_k,4.0);
  //bottom circle
  static CircleDomain cr_bo(vec_nk,vec_nk,8.0);

  //++ 2D domains ++

  static PlanarDomain sf_t(vec_k,vec_k), sf_b(vec_nk,vec_nk);
  static CylinderDomain cy(8.0,vec_k,vec_k);
  static SphericalDomain sp(vec_c,5.0);

  MeshDomain* base_domains[] = {
    &cr_to, &cr_ti, &cr_bo,
    &sf_t, &sf_b, &cy, &sp
  };
  const int NDOM = sizeof(base_domains)/sizeof(base_domains[0]);

  int dim_array[NDOM] = {
    1, 1, 1, 
    2, 2, 2, 2
  };

  //++ MESH & DOMAIN ++

  mesh.read_vtk(filename, err); MSQ_ERRZERO(err);
  DomainClassifier::classify_skin_geometrically (geom, &mesh, 0.001, base_domains, dim_array, NDOM, err);
  MSQ_ERRZERO(err);
  mesh.set_skin_flags( false, false, true, err ); MSQ_ERRZERO(err);
  
  return filename;
}

void create_input_mesh( const MeshParams& p, 
                        bool all_fixed,
                        MeshImpl& mesh, MsqError& err )
{
  const double z = 0;
  const int F = all_fixed;
  std::ofstream vtkfile( temp_file );
  double bx = all_fixed ? 0.5 * p.w : p.x;
  double by = all_fixed ? 0.5 * p.h : p.y;
  vtkfile << "# vtk DataFile Version 3.0" << std::endl
          << "Mesquite High Aspect Ratio test" << std::endl
          << "ASCII" << std::endl
          << "DATASET UNSTRUCTURED_GRID" << std::endl
          << "POINTS 9 float" << std::endl
          << 0.0 << ' ' << 0.0 << ' ' << z << std::endl
          << bx  << ' ' << 0.0 << ' ' << z << std::endl
          << p.w << ' ' << 0.0 << ' ' << z << std::endl
          << 0.0 << ' ' << by  << ' ' << z << std::endl
          << p.x << ' ' << p.y << ' ' << z << std::endl
          << p.w << ' ' << by  << ' ' << z << std::endl
          << 0.0 << ' ' << p.h << ' ' << z << std::endl
          << bx  << ' ' << p.h << ' ' << z << std::endl
          << p.w << ' ' << p.h << ' ' << z << std::endl
          << "CELLS 4 20" << std::endl
          << "4 0 1 4 3" << std::endl
          << "4 1 2 5 4" << std::endl
          << "4 4 5 8 7" << std::endl
          << "4 3 4 7 6" << std::endl
          << "CELL_TYPES 4" << std::endl
          << "9 9 9 9" << std::endl
          << "POINT_DATA 9" << std::endl
          << "SCALARS fixed int" << std::endl
          << "LOOKUP_TABLE default" << std::endl
          << "1 " << F << " 1" << std::endl
          <<  F << " 0 " << F  << std::endl
          << "1 " << F << " 1" << std::endl
          ;
  vtkfile.close();
  mesh.read_vtk( temp_file, err );
  std::remove( temp_file );
  MSQ_CHKERR(err);
}

  void VtkTest::test_read_unstructured( const char* filename )
  {
    MeshImpl mesh;
    MsqPrintError err(cout);
    
    mesh.read_vtk( filename, err );
    ASSERT_NO_ERROR(err);
    
      // Get mesh data
    std::vector<Mesh::VertexHandle> conn;
    std::vector<Mesh::ElementHandle> elems(38);
    std::vector<size_t> offsets(39);
    mesh.get_all_elements( elems, err );
    ASSERT_NO_ERROR(err);
    CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)38 );
    mesh.elements_get_attached_vertices( arrptr(elems), elems.size(),
                                         conn, offsets, err );
    ASSERT_NO_ERROR(err);

    unsigned i;
    struct meshdata { EntityTopology type; size_t nodes; size_t count; };
    meshdata list[] = {
      { Mesquite::HEXAHEDRON,    8,  8 },
      { Mesquite::QUADRILATERAL, 4,  4 },
      { Mesquite::PYRAMID,       5,  4 },
      { Mesquite::TETRAHEDRON,   4,  6 },
      { Mesquite::TRIANGLE,      3, 14 },
      { Mesquite::PRISM,         6,  2 }, 
      { Mesquite::MIXED,         0,  0 } };
      
      // Count expected lenght of connectivity list
    size_t conn_len = 0;
    for (i = 0; list[i].nodes; ++i)
      conn_len += list[i].nodes * list[i].count;
    CPPUNIT_ASSERT_EQUAL( conn_len, conn.size() );
    
    check_8hex_block( mesh, conn.begin() );
    check_4quad_block( mesh, conn.begin() + 64 );
  }

void create_input_mesh( double mid_x , MeshImpl& mesh, MsqError& err )
{
  std::ofstream vtkfile( temp_file );
  vtkfile << "# vtk DataFile Version 3.0" << std::endl
          << "Mesquite Syncronous Boundary test" << std::endl
          << "ASCII" << std::endl
          << "DATASET UNSTRUCTURED_GRID" << std::endl
          << "POINTS 9 float" << std::endl
          << min_x << ' ' << max_y << ' ' << z << std::endl
          << mid_x << ' ' << max_y << ' ' << z << std::endl
          << max_x << ' ' << max_y << ' ' << z << std::endl
          << min_x << ' ' << mid_y << ' ' << z << std::endl
          << mid_x << ' ' << mid_y << ' ' << z << std::endl
          << max_x << ' ' << mid_y << ' ' << z << std::endl
          << min_x << ' ' << min_y << ' ' << z << std::endl
          << mid_x << ' ' << min_y << ' ' << z << std::endl
          << max_x << ' ' << min_y << ' ' << z << std::endl
          << "CELLS 4 20" << std::endl
          << "4 1 0 3 4" << std::endl
          << "4 2 1 4 5" << std::endl
          << "4 4 3 6 7" << std::endl
          << "4 5 4 7 8" << std::endl
          << "CELL_TYPES 4" << std::endl
          << "9 9 9 9" << std::endl
          << "POINT_DATA 9" << std::endl
          << "SCALARS fixed int" << std::endl
          << "LOOKUP_TABLE default" << std::endl
          << "1 0 1" << std::endl
          << "1 0 1" << std::endl
          << "1 0 1" << std::endl
          ;
  vtkfile.close();
           
  mesh.read_vtk( temp_file, err );
  remove( temp_file );
}