MOAB: Mesh Oriented datABase
(version 5.2.1)
|
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