MOAB: Mesh Oriented datABase  (version 5.4.0)
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,
00407                             size_t num_elements,
00408                             const MBMesquite::MsqVertex* vtx_array,
00409                             size_t num_vertices );
00410 
00411     int tet_validity_check( const MBMesquite::MsqMeshEntity* element_array,
00412                             size_t num_elements,
00413                             const MBMesquite::MsqVertex* vtx_array );
00414 };
00415 // Check if the 2x2x2 brick of structured mesh
00416 // read from file is as expected.
00417 void VtkTest::check_8hex_structured( Mesh& mesh )
00418 {
00419     MsqPrintError err( cout );
00420 
00421     std::vector< Mesh::ElementHandle > elems( 8 );
00422     std::vector< Mesh::VertexHandle > verts( 64 );
00423     std::vector< size_t > offsets( 9 );
00424 
00425     mesh.get_all_elements( elems, err );
00426     ASSERT_NO_ERROR( err );
00427     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00428 
00429     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00430     ASSERT_NO_ERROR( err );
00431     CPPUNIT_ASSERT( verts.size() == 64 );
00432     CPPUNIT_ASSERT( offsets.size() == 9 );
00433 
00434     check_8hex_block( mesh, verts.begin() );
00435 }
00436 
00437 // Check if the 2x2x2 brick of hexes
00438 // read from file is as expected.
00439 void VtkTest::check_8hex_block( Mesh& mesh, std::vector< Mesh::VertexHandle >::iterator connectivity )
00440 {
00441     MsqPrintError err( cout );
00442     const int base_corners[]   = { 0, 1, 3, 4, 9, 10, 12, 13 };
00443     const int corner_offsets[] = { 0, 1, 4, 3, 9, 10, 13, 12 };
00444 
00445     for( int hex = 0; hex < 8; ++hex )
00446     {
00447         for( int node = 0; node < 8; ++node )
00448         {
00449             const int index = base_corners[hex] + corner_offsets[node];
00450             const int x     = index % 3;
00451             const int y     = ( index / 3 ) % 3;
00452             const int z     = index / 9;
00453             const Vector3D expected_coords( 1.5 * x, 1.5 * y, 1.5 * z );
00454             MsqVertex actual_coords;
00455             Mesh::VertexHandle* conn_ptr = &*connectivity;
00456             ++connectivity;
00457             mesh.vertices_get_coordinates( conn_ptr, &actual_coords, 1, err );
00458             CPPUNIT_ASSERT( !err );
00459             CPPUNIT_ASSERT( expected_coords.within_tolerance_box( actual_coords, DBL_EPSILON ) );
00460         }
00461     }
00462 }
00463 
00464 // Check if the 2x2 brick of structured mesh
00465 // read from file is as expected.
00466 void VtkTest::check_4quad_structured( Mesh& mesh )
00467 {
00468     MsqPrintError err( cout );
00469 
00470     std::vector< Mesh::ElementHandle > elems( 4 );
00471     std::vector< Mesh::VertexHandle > verts( 9 );
00472     std::vector< size_t > offsets( 5 );
00473 
00474     mesh.get_all_elements( elems, err );
00475     ASSERT_NO_ERROR( err );
00476     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)4 );
00477 
00478     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00479     ASSERT_NO_ERROR( err );
00480     CPPUNIT_ASSERT( verts.size() == 16 );
00481     CPPUNIT_ASSERT( offsets.size() == 5 );
00482 
00483     check_4quad_block( mesh, verts.begin() );
00484 }
00485 
00486 // Check if the 2x2 brick of quads
00487 // read from file is as expected.
00488 void VtkTest::check_4quad_block( Mesh& mesh, std::vector< Mesh::VertexHandle >::iterator connectivity )
00489 {
00490     MsqPrintError err( cout );
00491     const int base_corners[]   = { 0, 1, 3, 4 };
00492     const int corner_offsets[] = { 0, 1, 4, 3 };
00493     for( int quad = 0; quad < 4; ++quad )
00494     {
00495         for( int node = 0; node < 4; ++node )
00496         {
00497             const int index = base_corners[quad] + corner_offsets[node];
00498             const int x     = index % 3;
00499             const int y     = index / 3;
00500             const Vector3D expected_coords( 1.5 * x, 1.5 * y, 0.0 );
00501             MsqVertex actual_coords;
00502             Mesh::VertexHandle* conn_ptr = &*connectivity;
00503             ++connectivity;
00504             mesh.vertices_get_coordinates( conn_ptr, &actual_coords, 1, err );
00505             CPPUNIT_ASSERT( !err );
00506             CPPUNIT_ASSERT( expected_coords.within_tolerance_box( actual_coords, DBL_EPSILON ) );
00507         }
00508     }
00509 }
00510 
00511 // Test reading VTK unstructured mesh
00512 void VtkTest::test_read_unstructured()
00513 {
00514     FILE* file = fopen( temp_file_name, "w+" );
00515     CPPUNIT_ASSERT( file );
00516     int rval = fputs( mixed_unstructured_data, file );
00517     fclose( file );
00518     if( rval == EOF ) remove( temp_file_name );
00519     CPPUNIT_ASSERT( rval != EOF );
00520 
00521     test_read_unstructured( temp_file_name );
00522 }
00523 
00524 void VtkTest::test_read_unstructured( const char* filename )
00525 {
00526     MeshImpl mesh;
00527     MsqPrintError err( cout );
00528 
00529     mesh.read_vtk( filename, err );
00530     ASSERT_NO_ERROR( err );
00531 
00532     // Get mesh data
00533     std::vector< Mesh::VertexHandle > conn;
00534     std::vector< Mesh::ElementHandle > elems( 38 );
00535     std::vector< size_t > offsets( 39 );
00536     mesh.get_all_elements( elems, err );
00537     ASSERT_NO_ERROR( err );
00538     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)38 );
00539     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), conn, offsets, err );
00540     ASSERT_NO_ERROR( err );
00541 
00542     unsigned i;
00543     struct meshdata
00544     {
00545         EntityTopology type;
00546         size_t nodes;
00547         size_t count;
00548     };
00549     meshdata list[] = { { MBMesquite::HEXAHEDRON, 8, 8 }, { MBMesquite::QUADRILATERAL, 4, 4 },
00550                         { MBMesquite::PYRAMID, 5, 4 },    { MBMesquite::TETRAHEDRON, 4, 6 },
00551                         { MBMesquite::TRIANGLE, 3, 14 },  { MBMesquite::PRISM, 6, 2 },
00552                         { MBMesquite::MIXED, 0, 0 } };
00553 
00554     // Count expected lenght of connectivity list
00555     size_t conn_len = 0;
00556     for( i = 0; list[i].nodes; ++i )
00557         conn_len += list[i].nodes * list[i].count;
00558     CPPUNIT_ASSERT_EQUAL( conn_len, conn.size() );
00559 
00560     check_8hex_block( mesh, conn.begin() );
00561     check_4quad_block( mesh, conn.begin() + 64 );
00562 }
00563 
00564 // Test reading 2D Vtk structured-points mesh
00565 void VtkTest::test_read_structured_2d_points()
00566 {
00567     MeshImpl mesh;
00568     MsqPrintError err( cout );
00569 
00570     FILE* file = fopen( temp_file_name, "w+" );
00571     fputs( structured_2d_points_data, file );
00572     fclose( file );
00573 
00574     mesh.read_vtk( temp_file_name, err );
00575     remove( temp_file_name );
00576     ASSERT_NO_ERROR( err );
00577 
00578     check_4quad_structured( mesh );
00579 }
00580 
00581 // Test reading 3D Vtk structured-points mesh
00582 void VtkTest::test_read_structured_3d_points()
00583 {
00584     MeshImpl mesh;
00585     MsqPrintError err( cout );
00586 
00587     FILE* file = fopen( temp_file_name, "w+" );
00588     fputs( structured_3d_points_data, file );
00589     fclose( file );
00590 
00591     mesh.read_vtk( temp_file_name, err );
00592     remove( temp_file_name );
00593     ASSERT_NO_ERROR( err );
00594 
00595     check_8hex_structured( mesh );
00596 }
00597 
00598 // Test reading 3D Vtk structured-grid mesh
00599 void VtkTest::test_read_structured_grid()
00600 {
00601     MeshImpl mesh;
00602     MsqPrintError err( cout );
00603 
00604     FILE* file = fopen( temp_file_name, "w+" );
00605     fputs( structured_grid_data, file );
00606     fclose( file );
00607 
00608     mesh.read_vtk( temp_file_name, err );
00609     remove( temp_file_name );
00610     ASSERT_NO_ERROR( err );
00611 
00612     check_8hex_structured( mesh );
00613 }
00614 
00615 // Test reading 3D Vtk rectilinear-grid mesh
00616 void VtkTest::test_read_rectilinear_grid()
00617 {
00618     MeshImpl mesh;
00619     MsqPrintError err( cout );
00620 
00621     FILE* file = fopen( temp_file_name, "w+" );
00622     fputs( rectilinear_grid_data, file );
00623     fclose( file );
00624 
00625     mesh.read_vtk( temp_file_name, err );
00626     remove( temp_file_name );
00627     ASSERT_NO_ERROR( err );
00628 
00629     check_8hex_structured( mesh );
00630 }
00631 
00632 // Test reading Vtk simple (one-component) scalar attribute
00633 void VtkTest::test_read_simple_scalar_attrib()
00634 {
00635     MeshImpl mesh;
00636     MsqPrintError err( cout );
00637 
00638     FILE* file = fopen( temp_file_name, "w+" );
00639     fputs( structured_3d_points_data, file );
00640     fputs( simple_scalar_attrib, file );
00641     fclose( file );
00642 
00643     mesh.read_vtk( temp_file_name, err );
00644     remove( temp_file_name );
00645     ASSERT_NO_ERROR( err );
00646 
00647     std::vector< Mesh::ElementHandle > elems;
00648     mesh.get_all_elements( elems, err );
00649     CPPUNIT_ASSERT( !err );
00650     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00651 
00652     void* th = mesh.tag_get( "global_id", err );
00653     CPPUNIT_ASSERT( !err );
00654 
00655     std::string name;
00656     Mesh::TagType type;
00657     unsigned tagsize;
00658     mesh.tag_properties( th, name, type, tagsize, err );
00659     CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
00660 
00661     int elem_data[8];
00662     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_data, err );
00663     CPPUNIT_ASSERT( !err );
00664 
00665     for( int i = 0; i < 8; ++i )
00666         CPPUNIT_ASSERT( elem_data[i] == ( 1 + i ) );
00667 }
00668 
00669 // Test reading Vtk vector attribute
00670 void VtkTest::test_read_vector_attrib()
00671 {
00672     MeshImpl mesh;
00673     MsqPrintError err( cout );
00674 
00675     FILE* file = fopen( temp_file_name, "w+" );
00676     fputs( structured_3d_points_data, file );
00677     fputs( simple_vector_attrib, file );
00678     fclose( file );
00679 
00680     mesh.read_vtk( temp_file_name, err );
00681     remove( temp_file_name );
00682     ASSERT_NO_ERROR( err );
00683 
00684     std::vector< Mesh::ElementHandle > elems;
00685     mesh.get_all_elements( elems, err );
00686     CPPUNIT_ASSERT( !err );
00687     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00688 
00689     void* th = mesh.tag_get( "hexvect", err );
00690     CPPUNIT_ASSERT( !err );
00691 
00692     std::string name;
00693     Mesh::TagType type;
00694     unsigned tagsize;
00695     mesh.tag_properties( th, name, type, tagsize, err );
00696     CPPUNIT_ASSERT( !err && type == Mesh::DOUBLE && tagsize == 3 );
00697 
00698     double elem_data[24];
00699     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_data, err );
00700     CPPUNIT_ASSERT( !err );
00701 
00702     for( int i = 0; i < 8; ++i )
00703         CPPUNIT_ASSERT( Vector3D( elem_data + 3 * i ) == Vector3D( i + 1, i + 1, i + 1 ) );
00704 }
00705 
00706 // Test reading Vtk field attribute
00707 void VtkTest::test_read_field_attrib()
00708 {
00709     FILE* file = fopen( temp_file_name, "w+" );
00710     fputs( structured_3d_points_data, file );
00711     fputs( simple_field_attrib, file );
00712     fclose( file );
00713     check_field_attrib( temp_file_name );
00714 }
00715 
00716 void VtkTest::check_field_attrib( const char* temp_file_name )
00717 {
00718     MeshImpl mesh;
00719     MsqPrintError err( cout );
00720 
00721     mesh.read_vtk( temp_file_name, err );
00722     remove( temp_file_name );
00723     ASSERT_NO_ERROR( err );
00724 
00725     std::vector< Mesh::ElementHandle > elems;
00726     mesh.get_all_elements( elems, err );
00727     CPPUNIT_ASSERT( !err );
00728     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00729 
00730     std::string name;
00731     Mesh::TagType type;
00732     unsigned tagsize;
00733 
00734     void* th = mesh.tag_get( "test_field elem_vects", err );
00735     CPPUNIT_ASSERT( !err );
00736 
00737     mesh.tag_properties( th, name, type, tagsize, err );
00738     CPPUNIT_ASSERT( !err && type == Mesh::DOUBLE && tagsize == 3 );
00739 
00740     double elem_data[24];
00741     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_data, err );
00742     CPPUNIT_ASSERT( !err );
00743 
00744     for( int i = 0; i < 8; ++i )
00745         CPPUNIT_ASSERT( Vector3D( elem_data + 3 * i ) == Vector3D( i + 1, i + 1, i + 1 ) );
00746 
00747     th = mesh.tag_get( "test_field elem_ids", err );
00748     CPPUNIT_ASSERT( !err );
00749 
00750     mesh.tag_properties( th, name, type, tagsize, err );
00751     CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
00752 
00753     int elem_ids[8];
00754     mesh.tag_get_element_data( th, 8, arrptr( elems ), elem_ids, err );
00755     CPPUNIT_ASSERT( !err );
00756 
00757     for( int i = 0; i < 8; ++i )
00758         CPPUNIT_ASSERT( elem_ids[i] == i + 1 );
00759 
00760     th = mesh.tag_get( "field1", err );
00761     CPPUNIT_ASSERT( !err );
00762 
00763     mesh.tag_properties( th, name, type, tagsize, err );
00764     CPPUNIT_ASSERT( !err && type == Mesh::INT && tagsize == 1 );
00765 
00766     int values[8];
00767     mesh.tag_get_element_data( th, 8, arrptr( elems ), values, err );
00768     CPPUNIT_ASSERT( !err );
00769 
00770     for( int i = 0; i < 8; ++i )
00771         CPPUNIT_ASSERT( values[i] == 8 - i );
00772 }
00773 
00774 // Test writing quadtratic elements
00775 void VtkTest::test_write_field_attrib()
00776 {
00777     MeshImpl mesh;
00778     MsqPrintError err( cout );
00779 
00780     // Create file containing unstructured mesh test case
00781     FILE* file = fopen( temp_file_name, "w+" );
00782     fputs( structured_3d_points_data, file );
00783     fputs( simple_field_attrib, file );
00784     fclose( file );
00785 
00786     // Read unstructured mesh file
00787     mesh.read_vtk( temp_file_name, err );
00788     remove( temp_file_name );
00789     ASSERT_NO_ERROR( err );
00790 
00791     // Write unstructured mesh file back out
00792     mesh.write_vtk( temp_file_name, err );
00793     if( err ) remove( temp_file_name );
00794     ASSERT_NO_ERROR( err );
00795 
00796     // Check if file contained expected mesh
00797     check_field_attrib( temp_file_name );
00798 }
00799 
00800 // Test reading MeshImpl boundary-vertex bit
00801 // from Vtk scalar attribute.
00802 void VtkTest::test_read_fixed_attrib()
00803 {
00804     MeshImpl mesh;
00805     MsqPrintError err( cout );
00806 
00807     FILE* file = fopen( temp_file_name, "w+" );
00808     fputs( structured_3d_points_data, file );
00809     fputs( fixed_vertex_attrib, file );
00810     fclose( file );
00811 
00812     mesh.read_vtk( temp_file_name, err );
00813     remove( temp_file_name );
00814     ASSERT_NO_ERROR( err );
00815 
00816     std::vector< Mesh::ElementHandle > elems;
00817     mesh.get_all_elements( elems, err );
00818     CPPUNIT_ASSERT( !err );
00819     CPPUNIT_ASSERT_EQUAL( elems.size(), (size_t)8 );
00820 
00821     std::vector< Mesh::VertexHandle > verts;
00822     std::vector< size_t > offsets;
00823     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), verts, offsets, err );
00824     ASSERT_NO_ERROR( err );
00825 
00826     // get unique list of vertices
00827     std::vector< Mesh::VertexHandle >::iterator new_end;
00828     std::sort( verts.begin(), verts.end() );
00829     new_end = std::unique( verts.begin(), verts.end() );
00830     verts.resize( new_end - verts.begin() );
00831     CPPUNIT_ASSERT_EQUAL( verts.size(), (size_t)27 );
00832 
00833     // get fixed flag
00834     std::vector< bool > fixed;
00835     mesh.vertices_get_fixed_flag( arrptr( verts ), fixed, verts.size(), err );
00836     ASSERT_NO_ERROR( err );
00837     CPPUNIT_ASSERT_EQUAL( verts.size(), fixed.size() );
00838 
00839     for( int i = 0; i < 27; ++i )
00840     {
00841         bool should_be_fixed = ( i != 13 );
00842         CPPUNIT_ASSERT_EQUAL( should_be_fixed, (bool)fixed[i] );
00843     }
00844 }
00845 
00846 // Test writing VTK unstructured mesh
00847 void VtkTest::test_write()
00848 {
00849     MeshImpl mesh;
00850     MsqPrintError err( cout );
00851 
00852     // Create file containing unstructured mesh test case
00853     FILE* file = fopen( temp_file_name, "w+" );
00854     CPPUNIT_ASSERT( file );
00855     int rval = fputs( mixed_unstructured_data, file );
00856     fclose( file );
00857     if( rval == EOF ) remove( temp_file_name );
00858     CPPUNIT_ASSERT( rval != EOF );
00859 
00860     // Read unstructured mesh file
00861     mesh.read_vtk( temp_file_name, err );
00862     remove( temp_file_name );
00863     ASSERT_NO_ERROR( err );
00864 
00865     // Write unstructured mesh file back out
00866     mesh.write_vtk( temp_file_name, err );
00867     if( err ) remove( temp_file_name );
00868     ASSERT_NO_ERROR( err );
00869 
00870     // Check if file contained expected mesh
00871     test_read_unstructured( temp_file_name );
00872     remove( temp_file_name );
00873 }
00874 
00875 // Test reading VTK unstructured mesh quadratic elements
00876 void VtkTest::test_read_quadratic()
00877 {
00878     FILE* file = fopen( temp_file_name, "w+" );
00879     CPPUNIT_ASSERT( file );
00880     int rval = fputs( quadratic_unstructured_data, file );
00881     fclose( file );
00882     if( rval == EOF ) remove( temp_file_name );
00883     CPPUNIT_ASSERT( rval != EOF );
00884 
00885     test_read_quadratic( temp_file_name );
00886     remove( temp_file_name );
00887 }
00888 
00889 void VtkTest::test_read_quadratic( const char* filename )
00890 {
00891     const size_t NUM_ELEM = 8;
00892 
00893     MeshImpl mesh;
00894     MsqPrintError err( cout );
00895 
00896     mesh.read_vtk( filename, err );
00897     ASSERT_NO_ERROR( err );
00898 
00899     std::vector< Mesh::ElementHandle > elems( NUM_ELEM );
00900     mesh.get_all_elements( elems, err );
00901     ASSERT_NO_ERROR( err );
00902     CPPUNIT_ASSERT_EQUAL( elems.size(), NUM_ELEM );
00903 
00904     std::vector< Mesh::VertexHandle > conn;
00905     std::vector< size_t > offsets;
00906     mesh.elements_get_attached_vertices( arrptr( elems ), elems.size(), conn, offsets, err );
00907     ASSERT_NO_ERROR( err );
00908     CPPUNIT_ASSERT_EQUAL( conn.size(), (size_t)108 );
00909 
00910     EntityTopology types[NUM_ELEM];
00911     mesh.elements_get_topologies( arrptr( elems ), types, NUM_ELEM, err );
00912     ASSERT_NO_ERROR( err );
00913 
00914     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,
00915                                           1.0, -1.0, 1.0,  1.0, 1.0, 1.0,  -1.0, 1.0, 1.0,  -1.0, -1.0, 1.0 };
00916     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 };
00917     static const double pyr_corners[] = { 1.0,  -1.0, -1.0, 1.0,  1.0, -1.0, -1.0, 1.0,
00918                                           -1.0, -1.0, -1.0, -1.0, 0.0, 0.0,  1.0 };
00919     static const double pri_corners[] = { -1.0, -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, -1.0,
00920                                           -1.0, -1.0, 1.0,  1.0, 1.0, 1.0,  -1.0, 1.0, 1.0 };
00921     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 };
00922     static const unsigned tet_edges[] = { 0, 1, 1, 2, 2, 0, 0, 3, 1, 3, 2, 3 };
00923     static const unsigned pri_edges[] = { 0, 1, 1, 2, 2, 0, 0, 3, 1, 4, 2, 5, 3, 4, 4, 5, 5, 3 };
00924     static const unsigned pyr_edges[] = { 0, 1, 1, 2, 2, 3, 3, 0, 0, 4, 1, 4, 2, 4, 3, 4 };
00925     static const unsigned hex_faces[] = { 4, 0, 1, 5, 4, 4, 1, 2, 6, 5, 4, 2, 3, 7, 6,
00926                                           4, 3, 0, 4, 7, 4, 3, 2, 1, 0, 4, 4, 5, 6, 7 };
00927     static const struct
00928     {
00929         EntityTopology topology;
00930         unsigned num_corners;
00931         unsigned num_edges;
00932         unsigned num_faces;   // if non-zero expect mid-face nodes
00933         unsigned num_region;  // if non-zero expect mid-region node
00934         const double* corners;
00935         const unsigned* edges;
00936         const unsigned* faces;
00937     } expected_elems[NUM_ELEM] = { { MBMesquite::HEXAHEDRON, 8, 12, 0, 0, hex_corners, hex_edges, hex_faces },
00938                                    { MBMesquite::HEXAHEDRON, 8, 12, 6, 1, hex_corners, hex_edges, hex_faces },
00939                                    { MBMesquite::TETRAHEDRON, 4, 6, 0, 0, tet_corners, tet_edges, 0 },
00940                                    { MBMesquite::QUADRILATERAL, 4, 4, 0, 0, hex_corners, hex_edges, 0 },
00941                                    { MBMesquite::QUADRILATERAL, 4, 4, 0, 1, hex_corners, hex_edges, 0 },
00942                                    { MBMesquite::TRIANGLE, 3, 3, 0, 0, tet_corners, tet_edges, 0 },
00943                                    { MBMesquite::PRISM, 6, 9, 0, 0, pri_corners, pri_edges, 0 },
00944                                    { MBMesquite::PYRAMID, 5, 8, 0, 0, pyr_corners, pyr_edges, 0 } };
00945 
00946     MsqVertex have;
00947     std::vector< Mesh::VertexHandle >::iterator v_it = conn.begin();
00948     for( unsigned i = 0; i < NUM_ELEM; ++i )
00949     {
00950         CPPUNIT_ASSERT_EQUAL( expected_elems[i].topology, types[i] );
00951 
00952         size_t vtx_start = offsets[i];
00953         size_t vtx_end   = offsets[i + 1];
00954         size_t conn_len  = expected_elems[i].num_corners + expected_elems[i].num_edges + expected_elems[i].num_faces +
00955                           expected_elems[i].num_region;
00956         CPPUNIT_ASSERT_EQUAL( conn_len, vtx_end - vtx_start );
00957 
00958         for( unsigned c = 0; c < expected_elems[i].num_corners; ++c, ++v_it )
00959         {
00960             Vector3D expected( expected_elems[i].corners + 3 * c );
00961             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
00962             ASSERT_NO_ERROR( err );
00963             expected -= have;
00964             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
00965         }
00966 
00967         for( unsigned m = 0; m < expected_elems[i].num_edges; ++m, ++v_it )
00968         {
00969             unsigned start_idx = expected_elems[i].edges[2 * m];
00970             unsigned end_idx   = expected_elems[i].edges[2 * m + 1];
00971             Vector3D start( expected_elems[i].corners + 3 * start_idx );
00972             Vector3D end( expected_elems[i].corners + 3 * end_idx );
00973             Vector3D expected = 0.5 * ( start + end );
00974 
00975             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
00976             ASSERT_NO_ERROR( err );
00977 
00978             expected -= have;
00979             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
00980         }
00981 
00982         const unsigned* f_it = expected_elems[i].faces;
00983         for( unsigned m = 0; m < expected_elems[i].num_faces; ++m, ++v_it )
00984         {
00985             Vector3D expected( 0, 0, 0 );
00986             const unsigned face_size = *f_it;
00987             ++f_it;
00988             CPPUNIT_ASSERT( face_size == 3u || face_size == 4u );
00989             for( unsigned f = 0; f < face_size; ++f, ++f_it )
00990                 expected += Vector3D( expected_elems[i].corners + 3 * *f_it );
00991             expected /= face_size;
00992 
00993             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
00994             ASSERT_NO_ERROR( err );
00995 
00996             expected -= have;
00997             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
00998         }
00999 
01000         if( expected_elems[i].num_region )
01001         {
01002             CPPUNIT_ASSERT_EQUAL( 1u, expected_elems[i].num_region );
01003 
01004             Vector3D expected( 0, 0, 0 );
01005             for( unsigned m = 0; m < expected_elems[i].num_corners; ++m )
01006                 expected += Vector3D( expected_elems[i].corners + 3 * m );
01007             expected /= expected_elems[i].num_corners;
01008 
01009             mesh.vertices_get_coordinates( &*v_it, &have, 1, err );
01010             ASSERT_NO_ERROR( err );
01011 
01012             expected -= have;
01013             CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, expected.length(), DBL_EPSILON );
01014 
01015             ++v_it;
01016         }
01017     }
01018 }
01019 
01020 // Test writing quadtratic elements
01021 void VtkTest::test_write_quadratic()
01022 {
01023     MeshImpl mesh;
01024     MsqPrintError err( cout );
01025 
01026     // Create file containing unstructured mesh test case
01027     FILE* file = fopen( temp_file_name, "w+" );
01028     CPPUNIT_ASSERT( file );
01029     int rval = fputs( quadratic_unstructured_data, file );
01030     fclose( file );
01031     if( rval == EOF ) remove( temp_file_name );
01032     CPPUNIT_ASSERT( rval != EOF );
01033 
01034     // Read unstructured mesh file
01035     mesh.read_vtk( temp_file_name, err );
01036     remove( temp_file_name );
01037     ASSERT_NO_ERROR( err );
01038 
01039     // Write unstructured mesh file back out
01040     mesh.write_vtk( temp_file_name, err );
01041     if( err ) remove( temp_file_name );
01042     ASSERT_NO_ERROR( err );
01043 
01044     // Check if file contained expected mesh
01045     test_read_quadratic( temp_file_name );
01046     remove( temp_file_name );
01047 }
01048 
01049 void VtkTest::test_elements()
01050 {
01051     MBMesquite::MsqPrintError err( cout );
01052     MeshImpl mMesh;
01053     mMesh.read_vtk( MESH_FILES_DIR "2D/vtk/tris/untangled/equil_tri2.vtk", err );
01054     ASSERT_NO_ERROR( err );
01055     MBMesquite::MeshDomainAssoc mesh_and_domain = MBMesquite::MeshDomainAssoc( &mMesh, 0 );
01056     MBMesquite::Instruction::initialize_vertex_byte( &mesh_and_domain, 0, err );
01057     ASSERT_NO_ERROR( err );
01058 
01059     // Retrieve a patch
01060     MBMesquite::PatchData pd;
01061     pd.set_mesh( &mMesh );
01062     VertexPatches patch_set;
01063     patch_set.set_mesh( &mMesh );
01064     PatchIterator patches( &patch_set );
01065     patches.get_next_patch( pd, err );
01066     ASSERT_NO_ERROR( err );
01067 
01068     int free_vtx = pd.num_free_vertices();
01069     //    std::cout << "nb of free vertices: " << free_vtx << std::endl;
01070     CPPUNIT_ASSERT( free_vtx == 1 );
01071 
01072     MBMesquite::MsqMeshEntity* element_array = pd.get_element_array( err );
01073     ASSERT_NO_ERROR( err );
01074     size_t num_elements = pd.num_elements();
01075     CPPUNIT_ASSERT( num_elements == 6 );
01076 
01077     const MBMesquite::MsqVertex* vtx_array = pd.get_vertex_array( err );
01078     ASSERT_NO_ERROR( err );
01079     size_t num_vertices = pd.num_nodes();
01080     CPPUNIT_ASSERT( num_vertices == 7 );
01081 
01082     CPPUNIT_ASSERT( tri_check_validity( element_array, num_elements, vtx_array, num_vertices ) == 1 );
01083 
01084     patches.get_next_patch( pd, err );
01085     ASSERT_NO_ERROR( err );
01086 
01087     element_array = pd.get_element_array( err );
01088     ASSERT_NO_ERROR( err );
01089     num_elements = pd.num_elements();
01090     CPPUNIT_ASSERT( num_elements == 6 );
01091 
01092     vtx_array = pd.get_vertex_array( err );
01093     ASSERT_NO_ERROR( err );
01094     num_vertices = pd.num_nodes();
01095     CPPUNIT_ASSERT( num_vertices == 7 );
01096 
01097     CPPUNIT_ASSERT( tri_check_validity( element_array, num_elements, vtx_array, num_vertices ) == 1 );
01098 }
01099 
01100 int VtkTest::tri_check_validity( const MBMesquite::MsqMeshEntity* element_array,
01101                                  size_t num_elements,
01102                                  const MBMesquite::MsqVertex* vtx_array,
01103                                  size_t num_vertices )
01104 {
01105 
01106     /* check that the simplicial mesh is still valid,
01107        based on right handedness. Returns a 1 or a 0 */
01108     int valid   = 1;
01109     double dEps = 1.e-13;
01110 
01111     double x1, x2, x3, y1, y2, y3;  // z1, z2, z3;
01112     std::vector< size_t > vertex_indices;
01113 
01114     for( size_t i = 0; i < num_elements; i++ )
01115     {
01116         element_array[i].get_vertex_indices( vertex_indices );
01117 
01118         x1 = vtx_array[vertex_indices[0]][0];
01119         y1 = vtx_array[vertex_indices[0]][1];
01120         x2 = vtx_array[vertex_indices[1]][0];
01121         y2 = vtx_array[vertex_indices[1]][1];
01122         x3 = vtx_array[vertex_indices[2]][0];
01123         y3 = vtx_array[vertex_indices[2]][1];
01124 
01125         double a = x2 * y3 - x3 * y2;
01126         double b = y2 - y3;
01127         double c = x3 - x2;
01128 
01129         double area = .5 * ( a + b * x1 + c * y1 );
01130         if( area < dEps )
01131         {
01132             //          printf("x1 y1 = %f %f\n",x1,y1);
01133             //          printf("x2 y3 = %f %f\n",x2,y2);
01134             //          printf("x3 y3 = %f %f\n",x3,y3);
01135             //          printf("area = %f\n",area);
01136             valid = 0;
01137         }
01138     }
01139 
01140     return ( valid );
01141 }
01142 
01143 int VtkTest::tet_validity_check( const MBMesquite::MsqMeshEntity* element_array,
01144                                  size_t num_elements,
01145                                  const MBMesquite::MsqVertex* vtx_array )
01146 {
01147     int valid   = 1;
01148     double dEps = 1.e-13;
01149     double x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4;
01150     std::vector< size_t > vertex_indices;
01151 
01152     for( size_t i = 0; i < num_elements; i++ )
01153     {
01154         element_array[i].get_vertex_indices( vertex_indices );
01155 
01156         x1 = vtx_array[vertex_indices[0]][0];
01157         y1 = vtx_array[vertex_indices[0]][1];
01158         z1 = vtx_array[vertex_indices[0]][2];
01159 
01160         x2 = vtx_array[vertex_indices[1]][0];
01161         y2 = vtx_array[vertex_indices[1]][1];
01162         z2 = vtx_array[vertex_indices[1]][2];
01163 
01164         x3 = vtx_array[vertex_indices[2]][0];
01165         y3 = vtx_array[vertex_indices[2]][1];
01166         z3 = vtx_array[vertex_indices[2]][2];
01167 
01168         x4 = vtx_array[vertex_indices[3]][0];
01169         y4 = vtx_array[vertex_indices[3]][1];
01170         z4 = vtx_array[vertex_indices[3]][2];
01171 
01172         double dDX2 = x2 - x1;
01173         double dDX3 = x3 - x1;
01174         double dDX4 = x4 - x1;
01175 
01176         double dDY2 = y2 - y1;
01177         double dDY3 = y3 - y1;
01178         double dDY4 = y4 - y1;
01179 
01180         double dDZ2 = z2 - z1;
01181         double dDZ3 = z3 - z1;
01182         double dDZ4 = z4 - z1;
01183 
01184         /* dDet is proportional to the cell volume */
01185         double dDet = dDX2 * dDY3 * dDZ4 + dDX3 * dDY4 * dDZ2 + dDX4 * dDY2 * dDZ3 - dDZ2 * dDY3 * dDX4 -
01186                       dDZ3 * dDY4 * dDX2 - dDZ4 * dDY2 * dDX3;
01187 
01188         /* Compute a length scale based on edge lengths. */
01189         double dScale = ( sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) + ( z1 - z2 ) * ( z1 - z2 ) ) +
01190                           sqrt( ( x1 - x3 ) * ( x1 - x3 ) + ( y1 - y3 ) * ( y1 - y3 ) + ( z1 - z3 ) * ( z1 - z3 ) ) +
01191                           sqrt( ( x1 - x4 ) * ( x1 - x4 ) + ( y1 - y4 ) * ( y1 - y4 ) + ( z1 - z4 ) * ( z1 - z4 ) ) +
01192                           sqrt( ( x2 - x3 ) * ( x2 - x3 ) + ( y2 - y3 ) * ( y2 - y3 ) + ( z2 - z3 ) * ( z2 - z3 ) ) +
01193                           sqrt( ( x2 - x4 ) * ( x2 - x4 ) + ( y2 - y4 ) * ( y2 - y4 ) + ( z2 - z4 ) * ( z2 - z4 ) ) +
01194                           sqrt( ( x3 - x4 ) * ( x3 - x4 ) + ( y3 - y4 ) * ( y3 - y4 ) + ( z3 - z4 ) * ( z3 - z4 ) ) ) /
01195                         6.;
01196 
01197         /* Use the length scale to get a better idea if the tet is flat or
01198            just really small. */
01199         if( fabs( dScale ) < dEps )
01200         {
01201             return ( valid = 0 );
01202         }
01203         else
01204         {
01205             dDet /= ( dScale * dScale * dScale );
01206         }
01207 
01208         if( dDet > dEps )
01209         {
01210             valid = 1;
01211         }
01212         else if( dDet < -dEps )
01213         {
01214             valid = -1;
01215         }
01216         else
01217         {
01218             valid = 0;
01219         }
01220     }  // end for i=1,numElements
01221 
01222     return ( valid );
01223 }
01224 
01225 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( VtkTest, "VtkTest" );
01226 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( VtkTest, "Unit" );
01227 
01228 #else /* ifndef DEBUG */
01229 
01230 int main()
01231 {
01232     FILE* file;
01233 
01234     file = fopen( "2d_structured_points.vtk", "w" );
01235     fputs( structured_2d_points_data, file );
01236     fclose( file );
01237 
01238     file = fopen( "3d_structured_points.vtk", "w" );
01239     fputs( structured_3d_points_data, file );
01240     fputs( fixed_vertex_attrib, file );
01241     fclose( file );
01242 
01243     file = fopen( "structured_grid.vtk", "w" );
01244     fputs( structured_grid_data, file );
01245     fputs( fixed_vertex_attrib, file );
01246     fclose( file );
01247 
01248     file = fopen( "rectilinear_grid.vtk", "w" );
01249     fputs( rectilinear_grid_data, file );
01250     fputs( fixed_vertex_attrib, file );
01251     fclose( file );
01252 
01253     file = fopen( "mixed_unstructured.vtk", "w" );
01254     fputs( mixed_unstructured_data, file );
01255     fclose( file );
01256 
01257     file = fopen( "quadratic.vtk", "w" );
01258     fputs( quadratic_unstructured_data, file );
01259     fclose( file );
01260 
01261     return 0;
01262 }
01263 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines