MOAB: Mesh Oriented datABase  (version 5.2.1)
VtkTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2004 Sandia Corporation and Argonne National
00005     Laboratory.  Under the terms of Contract DE-AC04-94AL85000
00006     with Sandia Corporation, the U.S. Government retains certain
00007     rights in this software.
00008 
00009     This library is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU Lesser General Public
00011     License as published by the Free Software Foundation; either
00012     version 2.1 of the License, or (at your option) any later version.
00013 
00014     This library is distributed in the hope that it will be useful,
00015     but WITHOUT ANY WARRANTY; without even the implied warranty of
00016     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00017     Lesser General Public License for more details.
00018 
00019     You should have received a copy of the GNU Lesser General Public License
00020     (lgpl.txt) along with this library; if not, write to the Free Software
00021     Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
00022 
00023     diachin2@llnl.gov, djmelan@sandia.gov, mbrewer@sandia.gov,
00024     pknupp@sandia.gov, tleurent@mcs.anl.gov, tmunson@mcs.anl.gov
00025 
00026   ***************************************************************** */
00027 // -*- Mode : c++; tab-width: 2; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 2
00028 // -*-
00029 //
00030 
00031 #include "meshfiles.h"
00032 
00033 #include <iostream>
00034 using std::cout;
00035 
00036 #ifndef DEBUG
00037 #include "Mesquite.hpp"
00038 #include "PatchData.hpp"
00039 #include "MeshImpl.hpp"
00040 #include "MeshInterface.hpp"
00041 #include "VertexPatches.hpp"
00042 #include "PatchIterator.hpp"
00043 #include "UnitUtil.hpp"
00044 #include "Instruction.hpp"
00045 using MBMesquite::arrptr;
00046 using MBMesquite::EntityTopology;
00047 using MBMesquite::Mesh;
00048 using MBMesquite::MeshImpl;
00049 using MBMesquite::MsqPrintError;
00050 using MBMesquite::MsqVertex;
00051 using MBMesquite::PatchIterator;
00052 using MBMesquite::Vector3D;
00053 using MBMesquite::VertexPatches;
00054 #else
00055 #include <stdio.h>
00056 #endif
00057 
00058 #include <algorithm>
00059 
00060 extern const char temp_file_name[] = "VtkTest.vtk";
00061 
00062 // VTK file for 2x2x2 block of hexes as structured-point
00063 extern const char structured_3d_points_data[] = "# vtk DataFile Version 2.0\n"
00064                                                 "Mesquite Mesh\n"
00065                                                 "ASCII\n"
00066                                                 "DATASET STRUCTURED_POINTS\n"
00067                                                 "DIMENSIONS 3 3 3\n"
00068                                                 "ORIGIN 0 0 0\n"
00069                                                 "SPACING 1.5 1.5 1.5\n";
00070 
00071 // VTK file for 2x2 block of quads as structured-point
00072 extern const char structured_2d_points_data[] = "# vtk DataFile Version 2.0\n"
00073                                                 "Mesquite Mesh\n"
00074                                                 "ASCII\n"
00075                                                 "DATASET STRUCTURED_POINTS\n"
00076                                                 "DIMENSIONS 3 3 1\n"
00077                                                 "ORIGIN 0 0 0\n"
00078                                                 "SPACING 1.5 1.5 0.0\n";
00079 
00080 // VTK file for 2x2x2 block of hexes as structured-grid
00081 extern const char structured_grid_data[] = "# vtk DataFile Version 2.0\n"
00082                                            "Mesquite Mesh\n"
00083                                            "ASCII\n"
00084                                            "DATASET STRUCTURED_GRID\n"
00085                                            "DIMENSIONS 3 3 3\n"
00086                                            "POINTS 27 float\n"
00087                                            "0.0 0.0 0.0\n1.5 0.0 0.0\n3.0 0.0 0.0\n"
00088                                            "0.0 1.5 0.0\n1.5 1.5 0.0\n3.0 1.5 0.0\n"
00089                                            "0.0 3.0 0.0\n1.5 3.0 0.0\n3.0 3.0 0.0\n"
00090                                            "0.0 0.0 1.5\n1.5 0.0 1.5\n3.0 0.0 1.5\n"
00091                                            "0.0 1.5 1.5\n1.5 1.5 1.5\n3.0 1.5 1.5\n"
00092                                            "0.0 3.0 1.5\n1.5 3.0 1.5\n3.0 3.0 1.5\n"
00093                                            "0.0 0.0 3.0\n1.5 0.0 3.0\n3.0 0.0 3.0\n"
00094                                            "0.0 1.5 3.0\n1.5 1.5 3.0\n3.0 1.5 3.0\n"
00095                                            "0.0 3.0 3.0\n1.5 3.0 3.0\n3.0 3.0 3.0\n";
00096 
00097 // VTK file for 2x2x2 block of hexes as rectilinear-grid
00098 extern const char rectilinear_grid_data[] = "# vtk DataFile Version 2.0\n"
00099                                             "Mesquite Mesh\n"
00100                                             "ASCII\n"
00101                                             "DATASET RECTILINEAR_GRID\n"
00102                                             "DIMENSIONS 3 3 3\n"
00103                                             "X_COORDINATES 3 float\n"
00104                                             "0.0 1.5 3.0\n"
00105                                             "Y_COORDINATES 3 float\n"
00106                                             "0.0 1.5 3.0\n"
00107                                             "Z_COORDINATES 3 float\n"
00108                                             "0.0 1.5 3.0\n";
00109 
00110 // VTK file containing mixed-element unstructured mesh
00111 // First 8 elems result in same 2x2x2 block of hexes as
00112 // structured cases above.  The next 4 elems are the same
00113 // as the quads from the 2D structured-point file above.
00114 extern const char mixed_unstructured_data[] =
00115     "# vtk DataFile Version 2.0\n"
00116     "Mesquite Mesh\n"
00117     "ASCII\n"
00118     "DATASET UNSTRUCTURED_GRID\n"
00119     "POINTS 35 float\n"
00120     "\n"  // points for an 2x2x2 brick of hexes (same geom/topo as above structured meshes)
00121     "0.0 0.0 0.0\n1.5 0.0 0.0\n3.0 0.0 0.0\n"
00122     "0.0 1.5 0.0\n1.5 1.5 0.0\n3.0 1.5 0.0\n"
00123     "0.0 3.0 0.0\n1.5 3.0 0.0\n3.0 3.0 0.0\n"
00124     "0.0 0.0 1.5\n1.5 0.0 1.5\n3.0 0.0 1.5\n"
00125     "0.0 1.5 1.5\n1.5 1.5 1.5\n3.0 1.5 1.5\n"
00126     "0.0 3.0 1.5\n1.5 3.0 1.5\n3.0 3.0 1.5\n"
00127     "0.0 0.0 3.0\n1.5 0.0 3.0\n3.0 0.0 3.0\n"
00128     "0.0 1.5 3.0\n1.5 1.5 3.0\n3.0 1.5 3.0\n"
00129     "0.0 3.0 3.0\n1.5 3.0 3.0\n3.0 3.0 3.0\n"
00130     "\n"  // more points on +x side of brick for pyramids and tets
00131     "4 0.75 0.75\n4 2.25 0.75\n4 0.75 2.25\n4 2.25 2.25\n"
00132     "\n"  // more points for two prisms/wedges
00133     "6 0.75 0.75\n6 2.25 0.75\n6 0.75 2.25\n6 2.25 2.25\n"
00134     "\n"
00135     "CELLS 38 216\n"
00136     "\n"  // 8 hexes in 2x2x2 block
00137     "8  0  1  4  3  9 10 13 12\n"
00138     "8  1  2  5  4 10 11 14 13\n"
00139     "8  3  4  7  6 12 13 16 15\n"
00140     "8  4  5  8  7 13 14 17 16\n"
00141     "8  9 10 13 12 18 19 22 21\n"
00142     "8 10 11 14 13 19 20 23 22\n"
00143     "8 12 13 16 15 21 22 25 24\n"
00144     "8 13 14 17 16 22 23 26 25\n"
00145     "\n"  // Quads on -z face of hex (inverted to match structured data)
00146     "4 0 1 4 3\n"
00147     "4 1 2 5 4\n"
00148     "4 3 4 7 6\n"
00149     "4 4 5 8 7\n"
00150     "\n"  // some pyramids on the +x side of the block
00151     "5  2  5 14 11 27\n"
00152     "5  5  8 17 14 28\n"
00153     "5 11 14 23 20 29\n"
00154     "5 14 17 26 23 30\n"
00155     "\n"  // Some tetrahedrons around the pyramids
00156     "4  5 14 27 28\n"
00157     "4 14 28 17 30\n"
00158     "4 29 14 23 30\n"
00159     "4 11 27 14 29\n"
00160     "4 27 29 30 14\n"
00161     "4 28 27 30 14\n"
00162     "\n"  // Triangles bounding the pyramid/tet region
00163     "3  2  5 27\n"
00164     "3  5 28 27\n"
00165     "3  5  8 28\n"
00166     "3  8 17 28\n"
00167     "3 17 30 28\n"
00168     "3 17 26 30\n"
00169     "3 26 23 30\n"
00170     "3 23 29 30\n"
00171     "3 23 20 29\n"
00172     "3 20 11 29\n"
00173     "3 11 27 29\n"
00174     "3  2 27 11\n"
00175     "3 27 28 30\n"
00176     "3 27 30 29\n"
00177     "\n"  // A couple wedges/prisms
00178     "6 27 30 28 31 34 32\n"
00179     "6 30 27 29 34 31 33\n"
00180     "\n"
00181     "CELL_TYPES 38\n"
00182     "12 12 12 12 12 12 12 12\n"  //  8 hexes
00183     " 9  9  9  9\n"              //  4 quads
00184     "14 14 14 14\n"              //  4 pyramids
00185     "10 10 10 10 10 10\n"        //  6 tets
00186     " 5  5  5  5  5  5  5 "
00187     " 5  5  5  5  5  5  5\n"  // 14 tri
00188     "13 13\n";                // 2 wedges
00189 
00190 extern const char quadratic_unstructured_data[] =
00191     "# vtk DataFile Version 2.0\n"
00192     "Mesquite Mesh\n"
00193     "ASCII\n"
00194     "DATASET UNSTRUCTURED_GRID\n"
00195     "POINTS 41 float\n"
00196     "\n"                 // points for a single quadtratic hex, in ExodusII order
00197     " 1.0 -1.0 -1.0 \n"  //  0  bottom corners
00198     " 1.0  1.0 -1.0 \n"  //  1
00199     "-1.0  1.0 -1.0 \n"  //  2
00200     "-1.0 -1.0 -1.0 \n"  //  3
00201     " 1.0 -1.0  1.0 \n"  //  4  top corners
00202     " 1.0  1.0  1.0 \n"  //  5
00203     "-1.0  1.0  1.0 \n"  //  6
00204     "-1.0 -1.0  1.0 \n"  //  7
00205     " 1.0  0.0 -1.0 \n"  //  8  bottom mid-nodes
00206     " 0.0  1.0 -1.0 \n"  //  9
00207     "-1.0  0.0 -1.0 \n"  // 10
00208     " 0.0 -1.0 -1.0 \n"  // 11
00209     " 1.0 -1.0  0.0 \n"  // 12  side mid-nodes
00210     " 1.0  1.0  0.0 \n"  // 13
00211     "-1.0  1.0  0.0 \n"  // 14
00212     "-1.0 -1.0  0.0 \n"  // 15
00213     " 1.0  0.0  1.0 \n"  // 16  top mid-nodes
00214     " 0.0  1.0  1.0 \n"  // 17
00215     "-1.0  0.0  1.0 \n"  // 18
00216     " 0.0 -1.0  1.0 \n"  // 19
00217     " 1.0  0.0  0.0 \n"  // 20  mid-nodes for side faces
00218     " 0.0  1.0  0.0 \n"  // 21
00219     "-1.0  0.0  0.0 \n"  // 22
00220     " 0.0 -1.0  0.0 \n"  // 23
00221     " 0.0  0.0 -1.0 \n"  // 24  mid-nodes for top/bottom faces
00222     " 0.0  0.0  1.0 \n"  // 25
00223     " 0.0  0.0  0.0 \n"  // 26  mid-region
00224     "\n"                 // points for a single quadtatic tet, in ExodusII order
00225     " 1.0 -1.0 -1.0 \n"  // 27  base triangle
00226     " 1.0  1.0 -1.0 \n"  // 28
00227     "-1.0  0.0 -1.0 \n"  // 29
00228     " 0.0  0.0  1.0 \n"  // 30  apex
00229     " 1.0  0.0 -1.0 \n"  // 31  base mid-nodes
00230     " 0.0  0.5 -1.0 \n"  // 32
00231     " 0.0 -0.5 -1.0 \n"  // 33
00232     " 0.5 -0.5  0.0 \n"  // 34  side mid-nodes
00233     " 0.5  0.5  0.0 \n"  // 35
00234     "-0.5  0.0  0.0 \n"  // 36
00235     "\n"                 // mid-edge nodes for apex edges of quadratic pyramid
00236     " 0.5 -0.5 0 \n"     // 37
00237     " 0.5  0.5 0 \n"     // 38
00238     "-0.5  0.5 0 \n"     // 39
00239     "-0.5 -0.5 0 \n"     // 40
00240     "\n"
00241     "CELLS 8 116\n"
00242     "20 0 1 2 3 4 5 6 7 8 9 10 11 16 17 18 19 12 13 14 15\n"                       // Hex20
00243     "27 0 1 2 3 4 5 6 7 8 9 10 11 16 17 18 19 12 13 14 15 23 21 20 22 24 25 26\n"  // Hex27
00244     "10 27 28 29 30 31 32 33 34 35 36\n"                                           // Tet10
00245     "8 0 1 2 3 8 9 10 11\n"                                                        // Quad8 (base of hex)
00246     "9 0 1 2 3 8 9 10 11 24\n"                                                     // Quad9 (base of hex)
00247     "6 27 28 29 31 32 33\n"                                                        // Tri6 (base of tet)
00248     "15 3 2 1 7 6 5 10 9 24 18 17 25 15 14 13\n"                                   // quadratic prism as half of hex
00249     "13 0 1 2 3 25 8 9 10 11 37 38 39 40\n"  // quadratic pyramid with same base as hex
00250     "\n"
00251     "CELL_TYPES 8\n"
00252     "25 29 24 23 28 22 26 27\n";
00253 
00254 // A simple scalar attribute specifying point and hex
00255 // IDs.  May be appended to any of the above 3D structured
00256 // mesh files.
00257 extern const char simple_scalar_attrib[] = "CELL_DATA 8\n"
00258                                            "SCALARS global_id int 1\n"
00259                                            "LOOKUP_TABLE default\n"
00260                                            "1 2 3 4 5 6 7 8\n"
00261                                            "\n"
00262                                            "POINT_DATA 27\n"
00263                                            "SCALARS global_id int\n"
00264                                            "LOOKUP_TABLE default\n"
00265                                            " 1  2  3  4  5  6  7  8  9\n"
00266                                            "10 11 12 13 14 15 16 17 18\n"
00267                                            "19 20 21 22 23 24 25 26 27\n";
00268 
00269 // A VTK vector attribute.  May be appended to any of the
00270 // above 3D structured mesh files.
00271 extern const char simple_vector_attrib[] = "CELL_DATA 8\n"
00272                                            "VECTORS hexvect float\n"
00273                                            "1 1 1\n"
00274                                            "2 2 2\n"
00275                                            "3 3 3\n"
00276                                            "4 4 4\n"
00277                                            "5 5 5\n"
00278                                            "6 6 6\n"
00279                                            "7 7 7\n"
00280                                            "8 8 8\n";
00281 
00282 extern const char simple_field_attrib[] = "CELL_DATA 8\n"
00283                                           "FIELD test_field 2\n"
00284                                           "elem_ids 1 8 int\n"
00285                                           "1 2 3 4 5 6 7 8\n"
00286                                           "elem_vects 3 8 float\n"
00287                                           "1 1 1\n"
00288                                           "2 2 2\n"
00289                                           "3 3 3\n"
00290                                           "4 4 4\n"
00291                                           "5 5 5\n"
00292                                           "6 6 6\n"
00293                                           "7 7 7\n"
00294                                           "8 8 8\n"
00295                                           "FIELD field1 1\n"
00296                                           "values 1 8 int\n"
00297                                           "8 7 6 5 4 3 2 1\n";
00298 
00299 // A scalar VTK attribute with the name and datatype
00300 // expected by MeshImpl for specifying boundary vertices.
00301 // May be appended to any of the above 3D structured
00302 // mesh files
00303 extern const char fixed_vertex_attrib[] = "POINT_DATA 27\n"
00304                                           "SCALARS fixed float\n"
00305                                           "LOOKUP_TABLE default\n"
00306                                           "1 1 1 1 1 1 1 1 1\n"
00307                                           "1 1 1 1 0 1 1 1 1\n"
00308                                           "1 1 1 1 1 1 1 1 1\n";
00309 
00310 #ifndef DEBUG
00311 
00312 class VtkTest : public CppUnit::TestFixture
00313 {
00314   private:
00315     CPPUNIT_TEST_SUITE( VtkTest );
00316 
00317     // Original test for old Vtk parser
00318     CPPUNIT_TEST( test_elements );
00319 
00320     // Additional tests for new Vtk parser - J.Kraftcheck, 2004-10-12
00321     CPPUNIT_TEST( test_read_unstructured );
00322     CPPUNIT_TEST( test_read_structured_2d_points );
00323     CPPUNIT_TEST( test_read_structured_3d_points );
00324     CPPUNIT_TEST( test_read_structured_grid );
00325     CPPUNIT_TEST( test_read_rectilinear_grid );
00326     CPPUNIT_TEST( test_read_simple_scalar_attrib );
00327     CPPUNIT_TEST( test_read_vector_attrib );
00328     CPPUNIT_TEST( test_read_field_attrib );
00329     CPPUNIT_TEST( test_write_field_attrib );
00330     CPPUNIT_TEST( test_read_fixed_attrib );
00331     CPPUNIT_TEST( test_read_quadratic );
00332     CPPUNIT_TEST( test_write_quadratic );
00333 
00334     // Test writer - J.Kraftcheck, 2004-10-12
00335     CPPUNIT_TEST( test_write );
00336 
00337     CPPUNIT_TEST_SUITE_END();
00338 
00339   public:
00340     /* Automatically called by CppUnit before each test function. */
00341     void setUp() {}
00342 
00343     // Automatically called by CppUnit after each test function.
00344     void tearDown() {}
00345 
00346   public:
00347     VtkTest() {}
00348 
00349     // Check if the 2x2x2 brick of structured mesh
00350     // read from file is as expected.
00351     void check_8hex_structured( Mesh& mesh );
00352 
00353     // Check if the 2x2x2 brick of hexes
00354     // read from file is as expected.
00355     void check_8hex_block( Mesh& mesh, std::vector< Mesh::VertexHandle >::iterator connectivity );
00356 
00357     // Check if the 2x2 brick of structured mesh
00358     // read from file is as expected.
00359     void check_4quad_structured( Mesh& mesh );
00360 
00361     // Check if the 2x2 brick of quads
00362     // read from file is as expected.
00363     void check_4quad_block( Mesh& mesh, std::vector< Mesh::VertexHandle >::iterator connectivity );
00364 
00365     // Test reading VTK unstructured mesh
00366     void test_read_unstructured();
00367 
00368     void test_read_unstructured( const char* filename );
00369 
00370     // Test reading 2D Vtk structured-points mesh
00371     void test_read_structured_2d_points();
00372 
00373     // Test reading 3D Vtk structured-points mesh
00374     void test_read_structured_3d_points();
00375 
00376     // Test reading 3D Vtk structured-grid mesh
00377     void test_read_structured_grid();
00378 
00379     // Test reading 3D Vtk rectilinear-grid mesh
00380     void test_read_rectilinear_grid();
00381 
00382     // Test reading Vtk simple (one-component) scalar attribute
00383     void test_read_simple_scalar_attrib();
00384 
00385     // Test reading Vtk vector attribute
00386     void test_read_vector_attrib();
00387 
00388     // Test reading VTK FIELD attribute data
00389     void test_read_field_attrib();
00390     void test_write_field_attrib();
00391     void check_field_attrib( const char* file );
00392 
00393     // Test reading MeshImpl boundary-vertex bit
00394     // from Vtk scalar attribute.
00395     void test_read_fixed_attrib();
00396 
00397     // Test writing VTK unstructured mesh
00398     void test_write();
00399 
00400     void test_read_quadratic();
00401     void test_read_quadratic( const char* filename );
00402     void test_write_quadratic();
00403 
00404     void test_elements();
00405 
00406     int tri_check_validity( const MBMesquite::MsqMeshEntity* element_array, size_t num_elements,
00407                             const MBMesquite::MsqVertex* vtx_array, size_t num_vertices );
00408 
00409     int tet_validity_check( const MBMesquite::MsqMeshEntity* element_array, size_t num_elements,
00410                             const MBMesquite::MsqVertex* vtx_array );
00411 };
00412 // Check if the 2x2x2 brick of structured mesh
00413 // read from file is as expected.
00414 void VtkTest::check_8hex_structured( Mesh& mesh )
00415 {
00416     MsqPrintError err( cout );
00417 
00418     std::vector< Mesh::ElementHandle > elems( 8 );
00419     std::vector< Mesh::VertexHandle > verts( 64 );
00420     std::vector< size_t > offsets( 9 );
00421 
00422     mesh.get_all_elements( elems, err );
00423     ASSERT_NO_ERROR( err );
00424     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00425 
00426     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00427     ASSERT_NO_ERROR( err );
00428     CPPUNIT_ASSERT( verts.size() == 64 );
00429     CPPUNIT_ASSERT( offsets.size() == 9 );
00430 
00431     check_8hex_block( mesh, verts.begin() );
00432 }
00433 
00434 // Check if the 2x2x2 brick of hexes
00435 // read from file is as expected.
00436 void VtkTest::check_8hex_block( Mesh& mesh, std::vector< Mesh::VertexHandle >::iterator connectivity )
00437 {
00438     MsqPrintError err( cout );
00439     const int base_corners[]   = { 0, 1, 3, 4, 9, 10, 12, 13 };
00440     const int corner_offsets[] = { 0, 1, 4, 3, 9, 10, 13, 12 };
00441 
00442     for( int hex = 0; hex < 8; ++hex )
00443     {
00444         for( int node = 0; node < 8; ++node )
00445         {
00446             const int index = base_corners[hex] + corner_offsets[node];
00447             const int x     = index % 3;
00448             const int y     = ( index / 3 ) % 3;
00449             const int z     = index / 9;
00450             const Vector3D expected_coords( 1.5 * x, 1.5 * y, 1.5 * z );
00451             MsqVertex actual_coords;
00452             Mesh::VertexHandle* conn_ptr = &*connectivity;
00453             ++connectivity;
00454             mesh.vertices_get_coordinates( conn_ptr, &actual_coords, 1, err );
00455             CPPUNIT_ASSERT( !err );
00456             CPPUNIT_ASSERT( expected_coords.within_tolerance_box( actual_coords, DBL_EPSILON ) );
00457         }
00458     }
00459 }
00460 
00461 // Check if the 2x2 brick of structured mesh
00462 // read from file is as expected.
00463 void VtkTest::check_4quad_structured( Mesh& mesh )
00464 {
00465     MsqPrintError err( cout );
00466 
00467     std::vector< Mesh::ElementHandle > elems( 4 );
00468     std::vector< Mesh::VertexHandle > verts( 9 );
00469     std::vector< size_t > offsets( 5 );
00470 
00471     mesh.get_all_elements( elems, err );
00472     ASSERT_NO_ERROR( err );
00473     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)4 );
00474 
00475     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00476     ASSERT_NO_ERROR( err );
00477     CPPUNIT_ASSERT( verts.size() == 16 );
00478     CPPUNIT_ASSERT( offsets.size() == 5 );
00479 
00480     check_4quad_block( mesh, verts.begin() );
00481 }
00482 
00483 // Check if the 2x2 brick of quads
00484 // read from file is as expected.
00485 void VtkTest::check_4quad_block( Mesh& mesh, std::vector< Mesh::VertexHandle >::iterator connectivity )
00486 {
00487     MsqPrintError err( cout );
00488     const int base_corners[]   = { 0, 1, 3, 4 };
00489     const int corner_offsets[] = { 0, 1, 4, 3 };
00490     for( int quad = 0; quad < 4; ++quad )
00491     {
00492         for( int node = 0; node < 4; ++node )
00493         {
00494             const int index = base_corners[quad] + corner_offsets[node];
00495             const int x     = index % 3;
00496             const int y     = index / 3;
00497             const Vector3D expected_coords( 1.5 * x, 1.5 * y, 0.0 );
00498             MsqVertex actual_coords;
00499             Mesh::VertexHandle* conn_ptr = &*connectivity;
00500             ++connectivity;
00501             mesh.vertices_get_coordinates( conn_ptr, &actual_coords, 1, err );
00502             CPPUNIT_ASSERT( !err );
00503             CPPUNIT_ASSERT( expected_coords.within_tolerance_box( actual_coords, DBL_EPSILON ) );
00504         }
00505     }
00506 }
00507 
00508 // Test reading VTK unstructured mesh
00509 void VtkTest::test_read_unstructured()
00510 {
00511     FILE* file = fopen( temp_file_name, "w+" );
00512     CPPUNIT_ASSERT( file );
00513     int rval = fputs( mixed_unstructured_data, file );
00514     fclose( file );
00515     if( rval == EOF ) remove( temp_file_name );
00516     CPPUNIT_ASSERT( rval != EOF );
00517 
00518     test_read_unstructured( temp_file_name );
00519 }
00520 
00521 void VtkTest::test_read_unstructured( const char* filename )
00522 {
00523     MeshImpl mesh;
00524     MsqPrintError err( cout );
00525 
00526     mesh.read_vtk( filename, err );
00527     ASSERT_NO_ERROR( err );
00528 
00529     // Get mesh data
00530     std::vector< Mesh::VertexHandle > conn;
00531     std::vector< Mesh::ElementHandle > elems( 38 );
00532     std::vector< size_t > offsets( 39 );
00533     mesh.get_all_elements( elems, err );
00534     ASSERT_NO_ERROR( err );
00535     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)38 );
00536     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), conn, offsets, err );
00537     ASSERT_NO_ERROR( err );
00538 
00539     unsigned i;
00540     struct meshdata
00541     {
00542         EntityTopology type;
00543         size_t nodes;
00544         size_t count;
00545     };
00546     meshdata list[] = { { MBMesquite::HEXAHEDRON, 8, 8 }, { MBMesquite::QUADRILATERAL, 4, 4 },
00547                         { MBMesquite::PYRAMID, 5, 4 },    { MBMesquite::TETRAHEDRON, 4, 6 },
00548                         { MBMesquite::TRIANGLE, 3, 14 },  { MBMesquite::PRISM, 6, 2 },
00549                         { MBMesquite::MIXED, 0, 0 } };
00550 
00551     // Count expected lenght of connectivity list
00552     size_t conn_len = 0;
00553     for( i = 0; list[i].nodes; ++i )
00554         conn_len += list[i].nodes * list[i].count;
00555     CPPUNIT_ASSERT_EQUAL( conn_len, conn.size() );
00556 
00557     check_8hex_block( mesh, conn.begin() );
00558     check_4quad_block( mesh, conn.begin() + 64 );
00559 }
00560 
00561 // Test reading 2D Vtk structured-points mesh
00562 void VtkTest::test_read_structured_2d_points()
00563 {
00564     MeshImpl mesh;
00565     MsqPrintError err( cout );
00566 
00567     FILE* file = fopen( temp_file_name, "w+" );
00568     fputs( structured_2d_points_data, file );
00569     fclose( file );
00570 
00571     mesh.read_vtk( temp_file_name, err );
00572     remove( temp_file_name );
00573     ASSERT_NO_ERROR( err );
00574 
00575     check_4quad_structured( mesh );
00576 }
00577 
00578 // Test reading 3D Vtk structured-points mesh
00579 void VtkTest::test_read_structured_3d_points()
00580 {
00581     MeshImpl mesh;
00582     MsqPrintError err( cout );
00583 
00584     FILE* file = fopen( temp_file_name, "w+" );
00585     fputs( structured_3d_points_data, file );
00586     fclose( file );
00587 
00588     mesh.read_vtk( temp_file_name, err );
00589     remove( temp_file_name );
00590     ASSERT_NO_ERROR( err );
00591 
00592     check_8hex_structured( mesh );
00593 }
00594 
00595 // Test reading 3D Vtk structured-grid mesh
00596 void VtkTest::test_read_structured_grid()
00597 {
00598     MeshImpl mesh;
00599     MsqPrintError err( cout );
00600 
00601     FILE* file = fopen( temp_file_name, "w+" );
00602     fputs( structured_grid_data, file );
00603     fclose( file );
00604 
00605     mesh.read_vtk( temp_file_name, err );
00606     remove( temp_file_name );
00607     ASSERT_NO_ERROR( err );
00608 
00609     check_8hex_structured( mesh );
00610 }
00611 
00612 // Test reading 3D Vtk rectilinear-grid mesh
00613 void VtkTest::test_read_rectilinear_grid()
00614 {
00615     MeshImpl mesh;
00616     MsqPrintError err( cout );
00617 
00618     FILE* file = fopen( temp_file_name, "w+" );
00619     fputs( rectilinear_grid_data, file );
00620     fclose( file );
00621 
00622     mesh.read_vtk( temp_file_name, err );
00623     remove( temp_file_name );
00624     ASSERT_NO_ERROR( err );
00625 
00626     check_8hex_structured( mesh );
00627 }
00628 
00629 // Test reading Vtk simple (one-component) scalar attribute
00630 void VtkTest::test_read_simple_scalar_attrib()
00631 {
00632     MeshImpl mesh;
00633     MsqPrintError err( cout );
00634 
00635     FILE* file = fopen( temp_file_name, "w+" );
00636     fputs( structured_3d_points_data, file );
00637     fputs( simple_scalar_attrib, file );
00638     fclose( file );
00639 
00640     mesh.read_vtk( temp_file_name, err );
00641     remove( temp_file_name );
00642     ASSERT_NO_ERROR( err );
00643 
00644     std::vector< Mesh::ElementHandle > elems;
00645     mesh.get_all_elements( elems, err );
00646     CPPUNIT_ASSERT( !err );
00647     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00648 
00649     void* th = mesh.tag_get( "global_id", err );
00650     CPPUNIT_ASSERT( !err );
00651 
00652     std::string name;
00653     Mesh::TagType type;
00654     unsigned tagsize;
00655     mesh.tag_properties( th, name, type, tagsize, err );
00656     CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
00657 
00658     int elem_data[8];
00659     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_data, err );
00660     CPPUNIT_ASSERT( !err );
00661 
00662     for( int i = 0; i < 8; ++i )
00663         CPPUNIT_ASSERT( elem_data[i] == ( 1 + i ) );
00664 }
00665 
00666 // Test reading Vtk vector attribute
00667 void VtkTest::test_read_vector_attrib()
00668 {
00669     MeshImpl mesh;
00670     MsqPrintError err( cout );
00671 
00672     FILE* file = fopen( temp_file_name, "w+" );
00673     fputs( structured_3d_points_data, file );
00674     fputs( simple_vector_attrib, file );
00675     fclose( file );
00676 
00677     mesh.read_vtk( temp_file_name, err );
00678     remove( temp_file_name );
00679     ASSERT_NO_ERROR( err );
00680 
00681     std::vector< Mesh::ElementHandle > elems;
00682     mesh.get_all_elements( elems, err );
00683     CPPUNIT_ASSERT( !err );
00684     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00685 
00686     void* th = mesh.tag_get( "hexvect", err );
00687     CPPUNIT_ASSERT( !err );
00688 
00689     std::string name;
00690     Mesh::TagType type;
00691     unsigned tagsize;
00692     mesh.tag_properties( th, name, type, tagsize, err );
00693     CPPUNIT_ASSERT( !err && type == Mesh::DOUBLE && tagsize == 3 );
00694 
00695     double elem_data[24];
00696     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_data, err );
00697     CPPUNIT_ASSERT( !err );
00698 
00699     for( int i = 0; i < 8; ++i )
00700         CPPUNIT_ASSERT( Vector3D( elem_data + 3 * i ) == Vector3D( i + 1, i + 1, i + 1 ) );
00701 }
00702 
00703 // Test reading Vtk field attribute
00704 void VtkTest::test_read_field_attrib()
00705 {
00706     FILE* file = fopen( temp_file_name, "w+" );
00707     fputs( structured_3d_points_data, file );
00708     fputs( simple_field_attrib, file );
00709     fclose( file );
00710     check_field_attrib( temp_file_name );
00711 }
00712 
00713 void VtkTest::check_field_attrib( const char* temp_file_name )
00714 {
00715     MeshImpl mesh;
00716     MsqPrintError err( cout );
00717 
00718     mesh.read_vtk( temp_file_name, err );
00719     remove( temp_file_name );
00720     ASSERT_NO_ERROR( err );
00721 
00722     std::vector< Mesh::ElementHandle > elems;
00723     mesh.get_all_elements( elems, err );
00724     CPPUNIT_ASSERT( !err );
00725     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00726 
00727     std::string name;
00728     Mesh::TagType type;
00729     unsigned tagsize;
00730 
00731     void* th = mesh.tag_get( "test_field elem_vects", err );
00732     CPPUNIT_ASSERT( !err );
00733 
00734     mesh.tag_properties( th, name, type, tagsize, err );
00735     CPPUNIT_ASSERT( !err && type == Mesh::DOUBLE && tagsize == 3 );
00736 
00737     double elem_data[24];
00738     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_data, err );
00739     CPPUNIT_ASSERT( !err );
00740 
00741     for( int i = 0; i < 8; ++i )
00742         CPPUNIT_ASSERT( Vector3D( elem_data + 3 * i ) == Vector3D( i + 1, i + 1, i + 1 ) );
00743 
00744     th = mesh.tag_get( "test_field elem_ids", err );
00745     CPPUNIT_ASSERT( !err );
00746 
00747     mesh.tag_properties( th, name, type, tagsize, err );
00748     CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
00749 
00750     int elem_ids[8];
00751     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_ids, err );
00752     CPPUNIT_ASSERT( !err );
00753 
00754     for( int i = 0; i < 8; ++i )
00755         CPPUNIT_ASSERT( elem_ids[i] == i + 1 );
00756 
00757     th = mesh.tag_get( "field1", err );
00758     CPPUNIT_ASSERT( !err );
00759 
00760     mesh.tag_properties( th, name, type, tagsize, err );
00761     CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
00762 
00763     int values[8];
00764     mesh.tag_get_element_data( th, 8, arrptr( elems ), values, err );
00765     CPPUNIT_ASSERT( !err );
00766 
00767     for( int i = 0; i < 8; ++i )
00768         CPPUNIT_ASSERT( values[i] == 8 - i );
00769 }
00770 
00771 // Test writing quadtratic elements
00772 void VtkTest::test_write_field_attrib()
00773 {
00774     MeshImpl mesh;
00775     MsqPrintError err( cout );
00776 
00777     // Create file containing unstructured mesh test case
00778     FILE* file = fopen( temp_file_name, "w+" );
00779     fputs( structured_3d_points_data, file );
00780     fputs( simple_field_attrib, file );
00781     fclose( file );
00782 
00783     // Read unstructured mesh file
00784     mesh.read_vtk( temp_file_name, err );
00785     remove( temp_file_name );
00786     ASSERT_NO_ERROR( err );
00787 
00788     // Write unstructured mesh file back out
00789     mesh.write_vtk( temp_file_name, err );
00790     if( err ) remove( temp_file_name );
00791     ASSERT_NO_ERROR( err );
00792 
00793     // Check if file contained expected mesh
00794     check_field_attrib( temp_file_name );
00795 }
00796 
00797 // Test reading MeshImpl boundary-vertex bit
00798 // from Vtk scalar attribute.
00799 void VtkTest::test_read_fixed_attrib()
00800 {
00801     MeshImpl mesh;
00802     MsqPrintError err( cout );
00803 
00804     FILE* file = fopen( temp_file_name, "w+" );
00805     fputs( structured_3d_points_data, file );
00806     fputs( fixed_vertex_attrib, file );
00807     fclose( file );
00808 
00809     mesh.read_vtk( temp_file_name, err );
00810     remove( temp_file_name );
00811     ASSERT_NO_ERROR( err );
00812 
00813     std::vector< Mesh::ElementHandle > elems;
00814     mesh.get_all_elements( elems, err );
00815     CPPUNIT_ASSERT( !err );
00816     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00817 
00818     std::vector< Mesh::VertexHandle > verts;
00819     std::vector< size_t > offsets;
00820     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00821     ASSERT_NO_ERROR( err );
00822 
00823     // get unique list of vertices
00824     std::vector< Mesh::VertexHandle >::iterator new_end;
00825     std::sort( verts.begin(), verts.end() );
00826     new_end = std::unique( verts.begin(), verts.end() );
00827     verts.resize( new_end - verts.begin() );
00828     CPPUNIT_ASSERT_EQUAL( verts.size(), (size_t)27 );
00829 
00830     // get fixed flag
00831     std::vector< bool > fixed;
00832     mesh.vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err );
00833     ASSERT_NO_ERROR( err );
00834     CPPUNIT_ASSERT_EQUAL( verts.size(), fixed.size() );
00835 
00836     for( int i = 0; i < 27; ++i )
00837     {
00838         bool should_be_fixed = ( i != 13 );
00839         CPPUNIT_ASSERT_EQUAL( should_be_fixed, (bool)fixed[i] );
00840     }
00841 }
00842 
00843 // Test writing VTK unstructured mesh
00844 void VtkTest::test_write()
00845 {
00846     MeshImpl mesh;
00847     MsqPrintError err( cout );
00848 
00849     // Create file containing unstructured mesh test case
00850     FILE* file = fopen( temp_file_name, "w+" );
00851     CPPUNIT_ASSERT( file );
00852     int rval = fputs( mixed_unstructured_data, file );
00853     fclose( file );
00854     if( rval == EOF ) remove( temp_file_name );
00855     CPPUNIT_ASSERT( rval != EOF );
00856 
00857     // Read unstructured mesh file
00858     mesh.read_vtk( temp_file_name, err );
00859     remove( temp_file_name );
00860     ASSERT_NO_ERROR( err );
00861 
00862     // Write unstructured mesh file back out
00863     mesh.write_vtk( temp_file_name, err );
00864     if( err ) remove( temp_file_name );
00865     ASSERT_NO_ERROR( err );
00866 
00867     // Check if file contained expected mesh
00868     test_read_unstructured( temp_file_name );
00869     remove( temp_file_name );
00870 }
00871 
00872 // Test reading VTK unstructured mesh quadratic elements
00873 void VtkTest::test_read_quadratic()
00874 {
00875     FILE* file = fopen( temp_file_name, "w+" );
00876     CPPUNIT_ASSERT( file );
00877     int rval = fputs( quadratic_unstructured_data, file );
00878     fclose( file );
00879     if( rval == EOF ) remove( temp_file_name );
00880     CPPUNIT_ASSERT( rval != EOF );
00881 
00882     test_read_quadratic( temp_file_name );
00883     remove( temp_file_name );
00884 }
00885 
00886 void VtkTest::test_read_quadratic( const char* filename )
00887 {
00888     const size_t NUM_ELEM = 8;
00889 
00890     MeshImpl mesh;
00891     MsqPrintError err( cout );
00892 
00893     mesh.read_vtk( filename, err );
00894     ASSERT_NO_ERROR( err );
00895 
00896     std::vector< Mesh::ElementHandle > elems( NUM_ELEM );
00897     mesh.get_all_elements( elems, err );
00898     ASSERT_NO_ERROR( err );
00899     CPPUNIT_ASSERT_EQUAL( elems.size(), NUM_ELEM );
00900 
00901     std::vector< Mesh::VertexHandle > conn;
00902     std::vector< size_t > offsets;
00903     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), conn, offsets, err );
00904     ASSERT_NO_ERROR( err );
00905     CPPUNIT_ASSERT_EQUAL( conn.size(), (size_t)108 );
00906 
00907     EntityTopology types[NUM_ELEM];
00908     mesh.elements_get_topologies( arrptr( elems ), types, NUM_ELEM, err );
00909     ASSERT_NO_ERROR( err );
00910 
00911     static const double hex_corners[] = { 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0, -1.0, -1.0, -1.0,
00912                                           1.0, -1.0, 1.0,  1.0, 1.0, 1.0,  -1.0, 1.0, 1.0,  -1.0, -1.0, 1.0 };
00913     static const double tet_corners[] = { 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 0.0, -1.0, 0.0, 0.0, 1.0 };
00914     static const double pyr_corners[] = { 1.0,  -1.0, -1.0, 1.0,  1.0, -1.0, -1.0, 1.0,
00915                                           -1.0, -1.0, -1.0, -1.0, 0.0, 0.0,  1.0 };
00916     static const double pri_corners[] = { -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0,
00917                                           -1.0, -1.0, 1.0,  1.0, 1.0, 1.0,  -1.0, 1.0, 1.0 };
00918     static const unsigned hex_edges[] = { 0, 1, 1, 2, 2, 3, 3, 0, 0, 4, 1, 5, 2, 6, 3, 7, 4, 5, 5, 6, 6, 7, 7, 4 };
00919     static const unsigned tet_edges[] = { 0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3 };
00920     static const unsigned pri_edges[] = { 0, 1, 1, 2, 2, 0, 0, 3, 1, 4, 2, 5, 3, 4, 4, 5, 5, 3 };
00921     static const unsigned pyr_edges[] = { 0, 1, 1, 2, 2, 3, 3, 0, 0, 4, 1, 4, 2, 4, 3, 4 };
00922     static const unsigned hex_faces[] = { 4, 0, 1, 5, 4, 4, 1, 2, 6, 5, 4, 2, 3, 7, 6,
00923                                           4, 3, 0, 4, 7, 4, 3, 2, 1, 0, 4, 4, 5, 6, 7 };
00924     static const struct
00925     {
00926         EntityTopology topology;
00927         unsigned num_corners;
00928         unsigned num_edges;
00929         unsigned num_faces;   // if non-zero expect mid-face nodes
00930         unsigned num_region;  // if non-zero expect mid-region node
00931         const double* corners;
00932         const unsigned* edges;
00933         const unsigned* faces;
00934     } expected_elems[NUM_ELEM] = { { MBMesquite::HEXAHEDRON, 8, 12, 0, 0, hex_corners, hex_edges, hex_faces },
00935                                    { MBMesquite::HEXAHEDRON, 8, 12, 6, 1, hex_corners, hex_edges, hex_faces },
00936                                    { MBMesquite::TETRAHEDRON, 4, 6, 0, 0, tet_corners, tet_edges, 0 },
00937                                    { MBMesquite::QUADRILATERAL, 4, 4, 0, 0, hex_corners, hex_edges, 0 },
00938                                    { MBMesquite::QUADRILATERAL, 4, 4, 0, 1, hex_corners, hex_edges, 0 },
00939                                    { MBMesquite::TRIANGLE, 3, 3, 0, 0, tet_corners, tet_edges, 0 },
00940                                    { MBMesquite::PRISM, 6, 9, 0, 0, pri_corners, pri_edges, 0 },
00941                                    { MBMesquite::PYRAMID, 5, 8, 0, 0, pyr_corners, pyr_edges, 0 } };
00942 
00943     MsqVertex have;
00944     std::vector< Mesh::VertexHandle >::iterator v_it = conn.begin();
00945     for( unsigned i = 0; i < NUM_ELEM; ++i )
00946     {
00947         CPPUNIT_ASSERT_EQUAL( expected_elems[i].topology, types[i] );
00948 
00949         size_t vtx_start = offsets[i];
00950         size_t vtx_end   = offsets[i + 1];
00951         size_t conn_len  = expected_elems[i].num_corners + expected_elems[i].num_edges + expected_elems[i].num_faces +
00952                           expected_elems[i].num_region;
00953         CPPUNIT_ASSERT_EQUAL( conn_len, vtx_end - vtx_start );
00954 
00955         for( unsigned c = 0; c < expected_elems[i].num_corners; ++c, ++v_it )
00956         {
00957             Vector3D expected( expected_elems[i].corners + 3 * c );
00958             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
00959             ASSERT_NO_ERROR( err );
00960             expected -= have;
00961             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
00962         }
00963 
00964         for( unsigned m = 0; m < expected_elems[i].num_edges; ++m, ++v_it )
00965         {
00966             unsigned start_idx = expected_elems[i].edges[2 * m];
00967             unsigned end_idx   = expected_elems[i].edges[2 * m + 1];
00968             Vector3D start( expected_elems[i].corners + 3 * start_idx );
00969             Vector3D end( expected_elems[i].corners + 3 * end_idx );
00970             Vector3D expected = 0.5 * ( start + end );
00971 
00972             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
00973             ASSERT_NO_ERROR( err );
00974 
00975             expected -= have;
00976             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
00977         }
00978 
00979         const unsigned* f_it = expected_elems[i].faces;
00980         for( unsigned m = 0; m < expected_elems[i].num_faces; ++m, ++v_it )
00981         {
00982             Vector3D expected( 0, 0, 0 );
00983             const unsigned face_size = *f_it;
00984             ++f_it;
00985             CPPUNIT_ASSERT( face_size == 3u || face_size == 4u );
00986             for( unsigned f = 0; f < face_size; ++f, ++f_it )
00987                 expected += Vector3D( expected_elems[i].corners + 3 * *f_it );
00988             expected /= face_size;
00989 
00990             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
00991             ASSERT_NO_ERROR( err );
00992 
00993             expected -= have;
00994             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
00995         }
00996 
00997         if( expected_elems[i].num_region )
00998         {
00999             CPPUNIT_ASSERT_EQUAL( 1u, expected_elems[i].num_region );
01000 
01001             Vector3D expected( 0, 0, 0 );
01002             for( unsigned m = 0; m < expected_elems[i].num_corners; ++m )
01003                 expected += Vector3D( expected_elems[i].corners + 3 * m );
01004             expected /= expected_elems[i].num_corners;
01005 
01006             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
01007             ASSERT_NO_ERROR( err );
01008 
01009             expected -= have;
01010             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
01011 
01012             ++v_it;
01013         }
01014     }
01015 }
01016 
01017 // Test writing quadtratic elements
01018 void VtkTest::test_write_quadratic()
01019 {
01020     MeshImpl mesh;
01021     MsqPrintError err( cout );
01022 
01023     // Create file containing unstructured mesh test case
01024     FILE* file = fopen( temp_file_name, "w+" );
01025     CPPUNIT_ASSERT( file );
01026     int rval = fputs( quadratic_unstructured_data, file );
01027     fclose( file );
01028     if( rval == EOF ) remove( temp_file_name );
01029     CPPUNIT_ASSERT( rval != EOF );
01030 
01031     // Read unstructured mesh file
01032     mesh.read_vtk( temp_file_name, err );
01033     remove( temp_file_name );
01034     ASSERT_NO_ERROR( err );
01035 
01036     // Write unstructured mesh file back out
01037     mesh.write_vtk( temp_file_name, err );
01038     if( err ) remove( temp_file_name );
01039     ASSERT_NO_ERROR( err );
01040 
01041     // Check if file contained expected mesh
01042     test_read_quadratic( temp_file_name );
01043     remove( temp_file_name );
01044 }
01045 
01046 void VtkTest::test_elements()
01047 {
01048     MBMesquite::MsqPrintError err( cout );
01049     MeshImpl mMesh;
01050     mMesh.read_vtk( MESH_FILES_DIR "2D/vtk/tris/untangled/equil_tri2.vtk", err );
01051     ASSERT_NO_ERROR( err );
01052     MBMesquite::MeshDomainAssoc mesh_and_domain = MBMesquite::MeshDomainAssoc( &mMesh, 0 );
01053     MBMesquite::Instruction::initialize_vertex_byte( &mesh_and_domain, 0, err );
01054     ASSERT_NO_ERROR( err );
01055 
01056     // Retrieve a patch
01057     MBMesquite::PatchData pd;
01058     pd.set_mesh( &mMesh );
01059     VertexPatches patch_set;
01060     patch_set.set_mesh( &mMesh );
01061     PatchIterator patches( &patch_set );
01062     patches.get_next_patch( pd, err );
01063     ASSERT_NO_ERROR( err );
01064 
01065     int free_vtx = pd.num_free_vertices();
01066     //    std::cout << "nb of free vertices: " << free_vtx << std::endl;
01067     CPPUNIT_ASSERT( free_vtx == 1 );
01068 
01069     MBMesquite::MsqMeshEntity* element_array = pd.get_element_array( err );
01070     ASSERT_NO_ERROR( err );
01071     size_t num_elements = pd.num_elements();
01072     CPPUNIT_ASSERT( num_elements == 6 );
01073 
01074     const MBMesquite::MsqVertex* vtx_array = pd.get_vertex_array( err );
01075     ASSERT_NO_ERROR( err );
01076     size_t num_vertices = pd.num_nodes();
01077     CPPUNIT_ASSERT( num_vertices == 7 );
01078 
01079     CPPUNIT_ASSERT( tri_check_validity( element_array, num_elements, vtx_array, num_vertices ) == 1 );
01080 
01081     patches.get_next_patch( pd, err );
01082     ASSERT_NO_ERROR( err );
01083 
01084     element_array = pd.get_element_array( err );
01085     ASSERT_NO_ERROR( err );
01086     num_elements = pd.num_elements();
01087     CPPUNIT_ASSERT( num_elements == 6 );
01088 
01089     vtx_array = pd.get_vertex_array( err );
01090     ASSERT_NO_ERROR( err );
01091     num_vertices = pd.num_nodes();
01092     CPPUNIT_ASSERT( num_vertices == 7 );
01093 
01094     CPPUNIT_ASSERT( tri_check_validity( element_array, num_elements, vtx_array, num_vertices ) == 1 );
01095 }
01096 
01097 int VtkTest::tri_check_validity( const MBMesquite::MsqMeshEntity* element_array, size_t num_elements,
01098                                  const MBMesquite::MsqVertex* vtx_array, size_t num_vertices )
01099 {
01100 
01101     /* check that the simplicial mesh is still valid,
01102        based on right handedness. Returns a 1 or a 0 */
01103     int valid   = 1;
01104     double dEps = 1.e-13;
01105 
01106     double x1, x2, x3, y1, y2, y3;  // z1, z2, z3;
01107     std::vector< size_t > vertex_indices;
01108 
01109     for( size_t i = 0; i < num_elements; i++ )
01110     {
01111         element_array[i].get_vertex_indices( vertex_indices );
01112 
01113         x1 = vtx_array[vertex_indices[0]][0];
01114         y1 = vtx_array[vertex_indices[0]][1];
01115         x2 = vtx_array[vertex_indices[1]][0];
01116         y2 = vtx_array[vertex_indices[1]][1];
01117         x3 = vtx_array[vertex_indices[2]][0];
01118         y3 = vtx_array[vertex_indices[2]][1];
01119 
01120         double a = x2 * y3 - x3 * y2;
01121         double b = y2 - y3;
01122         double c = x3 - x2;
01123 
01124         double area = .5 * ( a + b * x1 + c * y1 );
01125         if( area < dEps )
01126         {
01127             //          printf("x1 y1 = %f %f\n",x1,y1);
01128             //          printf("x2 y3 = %f %f\n",x2,y2);
01129             //          printf("x3 y3 = %f %f\n",x3,y3);
01130             //          printf("area = %f\n",area);
01131             valid = 0;
01132         }
01133     }
01134 
01135     return ( valid );
01136 }
01137 
01138 int VtkTest::tet_validity_check( const MBMesquite::MsqMeshEntity* element_array, size_t num_elements,
01139                                  const MBMesquite::MsqVertex* vtx_array )
01140 {
01141     int valid   = 1;
01142     double dEps = 1.e-13;
01143     double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
01144     std::vector< size_t > vertex_indices;
01145 
01146     for( size_t i = 0; i < num_elements; i++ )
01147     {
01148         element_array[i].get_vertex_indices( vertex_indices );
01149 
01150         x1 = vtx_array[vertex_indices[0]][0];
01151         y1 = vtx_array[vertex_indices[0]][1];
01152         z1 = vtx_array[vertex_indices[0]][2];
01153 
01154         x2 = vtx_array[vertex_indices[1]][0];
01155         y2 = vtx_array[vertex_indices[1]][1];
01156         z2 = vtx_array[vertex_indices[1]][2];
01157 
01158         x3 = vtx_array[vertex_indices[2]][0];
01159         y3 = vtx_array[vertex_indices[2]][1];
01160         z3 = vtx_array[vertex_indices[2]][2];
01161 
01162         x4 = vtx_array[vertex_indices[3]][0];
01163         y4 = vtx_array[vertex_indices[3]][1];
01164         z4 = vtx_array[vertex_indices[3]][2];
01165 
01166         double dDX2 = x2 - x1;
01167         double dDX3 = x3 - x1;
01168         double dDX4 = x4 - x1;
01169 
01170         double dDY2 = y2 - y1;
01171         double dDY3 = y3 - y1;
01172         double dDY4 = y4 - y1;
01173 
01174         double dDZ2 = z2 - z1;
01175         double dDZ3 = z3 - z1;
01176         double dDZ4 = z4 - z1;
01177 
01178         /* dDet is proportional to the cell volume */
01179         double dDet = dDX2 * dDY3 * dDZ4 + dDX3 * dDY4 * dDZ2 + dDX4 * dDY2 * dDZ3 - dDZ2 * dDY3 * dDX4 -
01180                       dDZ3 * dDY4 * dDX2 - dDZ4 * dDY2 * dDX3;
01181 
01182         /* Compute a length scale based on edge lengths. */
01183         double dScale = ( sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) + ( z1 - z2 ) * ( z1 - z2 ) ) +
01184                           sqrt( ( x1 - x3 ) * ( x1 - x3 ) + ( y1 - y3 ) * ( y1 - y3 ) + ( z1 - z3 ) * ( z1 - z3 ) ) +
01185                           sqrt( ( x1 - x4 ) * ( x1 - x4 ) + ( y1 - y4 ) * ( y1 - y4 ) + ( z1 - z4 ) * ( z1 - z4 ) ) +
01186                           sqrt( ( x2 - x3 ) * ( x2 - x3 ) + ( y2 - y3 ) * ( y2 - y3 ) + ( z2 - z3 ) * ( z2 - z3 ) ) +
01187                           sqrt( ( x2 - x4 ) * ( x2 - x4 ) + ( y2 - y4 ) * ( y2 - y4 ) + ( z2 - z4 ) * ( z2 - z4 ) ) +
01188                           sqrt( ( x3 - x4 ) * ( x3 - x4 ) + ( y3 - y4 ) * ( y3 - y4 ) + ( z3 - z4 ) * ( z3 - z4 ) ) ) /
01189                         6.;
01190 
01191         /* Use the length scale to get a better idea if the tet is flat or
01192            just really small. */
01193         if( fabs( dScale ) < dEps ) { return ( valid = 0 ); }
01194         else
01195         {
01196             dDet /= ( dScale * dScale * dScale );
01197         }
01198 
01199         if( dDet > dEps ) { valid = 1; }
01200         else if( dDet < -dEps )
01201         {
01202             valid = -1;
01203         }
01204         else
01205         {
01206             valid = 0;
01207         }
01208     }  // end for i=1,numElements
01209 
01210     return ( valid );
01211 }
01212 
01213 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( VtkTest, "VtkTest" );
01214 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( VtkTest, "Unit" );
01215 
01216 #else /* ifndef DEBUG */
01217 
01218 int main()
01219 {
01220     FILE* file;
01221 
01222     file = fopen( "2d_structured_points.vtk", "w" );
01223     fputs( structured_2d_points_data, file );
01224     fclose( file );
01225 
01226     file = fopen( "3d_structured_points.vtk", "w" );
01227     fputs( structured_3d_points_data, file );
01228     fputs( fixed_vertex_attrib, file );
01229     fclose( file );
01230 
01231     file = fopen( "structured_grid.vtk", "w" );
01232     fputs( structured_grid_data, file );
01233     fputs( fixed_vertex_attrib, file );
01234     fclose( file );
01235 
01236     file = fopen( "rectilinear_grid.vtk", "w" );
01237     fputs( rectilinear_grid_data, file );
01238     fputs( fixed_vertex_attrib, file );
01239     fclose( file );
01240 
01241     file = fopen( "mixed_unstructured.vtk", "w" );
01242     fputs( mixed_unstructured_data, file );
01243     fclose( file );
01244 
01245     file = fopen( "quadratic.vtk", "w" );
01246     fputs( quadratic_unstructured_data, file );
01247     fclose( file );
01248 
01249     return 0;
01250 }
01251 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines