MOAB: Mesh Oriented datABase  (version 5.2.1)
BoundedCylinderDomainTest.cpp
Go to the documentation of this file.
00001 #include "BoundedCylinderDomain.hpp"
00002 #include "MsqError.hpp"
00003 #include "MeshImpl.hpp"
00004 #include "MsqVertex.hpp"
00005 #include <cppunit/extensions/HelperMacros.h>
00006 #include <vector>
00007 #include <algorithm>
00008 #include <stdlib.h>
00009 #include <stdio.h>
00010 #include <iostream>
00011 
00012 const double EPSILON = 1e-6;
00013 #define ASSERT_VECTORS_EQUAL( A, B ) CPPUNIT_ASSERT( ( ( A ) - ( B ) ).length() < EPSILON )
00014 
00015 using namespace MBMesquite;
00016 
00017 class BoundedCylinderDomainTest : public CppUnit::TestFixture
00018 {
00019   private:
00020     CPPUNIT_TEST_SUITE( BoundedCylinderDomainTest );
00021     CPPUNIT_TEST( test_domain_DoF );
00022     CPPUNIT_TEST( test_z_snap_to );
00023     CPPUNIT_TEST( test_x_snap_to );
00024     CPPUNIT_TEST( test_create_curve_from_mesh );
00025     CPPUNIT_TEST_SUITE_END();
00026 
00027   public:
00028     BoundedCylinderDomainTest() {}
00029 
00030     void setUp() {}
00031     void tearDown() {}
00032 
00033     void test_domain_DoF();
00034 
00035     void test_z_snap_to();
00036     void test_x_snap_to();
00037 
00038     void test_create_curve_from_mesh();
00039 };
00040 
00041 void BoundedCylinderDomainTest::test_domain_DoF()
00042 {
00043     MsqError err;
00044     const Mesh::EntityHandle list[] = { (void*)5, (void*)9, (void*)1, (void*)2, (void*)3 };
00045     const size_t len                = sizeof( list ) / sizeof( list[0] );
00046 
00047     BoundedCylinderDomain d( 1000 );
00048     std::vector< Mesh::EntityHandle > handles;
00049     std::copy( list, list + len, std::back_inserter( handles ) );
00050     d.create_curve( 10, handles );
00051 
00052     const Mesh::EntityHandle list2[] = { (void*)1, (void*)2, (void*)3, (void*)4, (void*)5,
00053                                          (void*)6, (void*)7, (void*)8, (void*)9, (void*)10 };
00054     const size_t len2                = sizeof( list2 ) / sizeof( list2[0] );
00055     unsigned short dof_list[len2];
00056     d.domain_DoF( list2, dof_list, len2, err );
00057     CPPUNIT_ASSERT( !err );
00058 
00059     for( size_t i = 0; i < len2; ++i )
00060     {
00061         Mesh::EntityHandle handle   = list2[i];
00062         const Mesh::EntityHandle* f = std::find( list, list + len, handle );
00063         unsigned short expected_dim = ( f < list + len ) ? 1 : 2;
00064         CPPUNIT_ASSERT_EQUAL( expected_dim, dof_list[i] );
00065     }
00066 }
00067 
00068 void BoundedCylinderDomainTest::test_z_snap_to()
00069 {
00070     const Mesh::EntityHandle HANDLE = (void*)20;
00071     const double ZPLANE             = 5;
00072 
00073     // Define a circle in the Z=5 plane, centered at (0,0,5)
00074     // with a radius of 10.
00075     std::vector< Mesh::EntityHandle > handles( 1 );
00076     handles[0] = HANDLE;
00077     BoundedCylinderDomain d( 10 );
00078     d.create_curve( ZPLANE, handles );
00079 
00080     Vector3D point;
00081 
00082     // try a point inside and below the circle
00083     point.set( 5, 5, 2 );
00084     d.snap_to( HANDLE, point );
00085     CPPUNIT_ASSERT_DOUBLES_EQUAL( ZPLANE, point.z(), EPSILON );
00086     CPPUNIT_ASSERT_DOUBLES_EQUAL( point.x(), point.y(), EPSILON );
00087     point[2] = 0.0;
00088     CPPUNIT_ASSERT_DOUBLES_EQUAL( d.radius(), point.length(), EPSILON );
00089 
00090     // try a point outside and above the circle
00091     point.set( 50, 50, 50 );
00092     d.snap_to( HANDLE, point );
00093     CPPUNIT_ASSERT_DOUBLES_EQUAL( ZPLANE, point.z(), EPSILON );
00094     CPPUNIT_ASSERT_DOUBLES_EQUAL( point.x(), point.y(), EPSILON );
00095     point[2] = 0.0;
00096     CPPUNIT_ASSERT_DOUBLES_EQUAL( d.radius(), point.length(), EPSILON );
00097 
00098     // try a point on the circle
00099     Vector3D on_circ( d.radius(), 0, ZPLANE );
00100     point = on_circ;
00101     d.snap_to( HANDLE, point );
00102     ASSERT_VECTORS_EQUAL( on_circ, point );
00103 
00104     // try a point at the center of the circle
00105     point.set( d.radius(), 0, ZPLANE );
00106     d.snap_to( HANDLE, point );
00107     CPPUNIT_ASSERT_DOUBLES_EQUAL( ZPLANE, point.z(), EPSILON );
00108     point[2] = 0.0;
00109     CPPUNIT_ASSERT_DOUBLES_EQUAL( d.radius(), point.length(), EPSILON );
00110 
00111     // try a different handle to make sure a handle not
00112     // on the the circle results in a point not on the circle
00113     double pplane = ZPLANE - 5;
00114     point.set( 0, 0, pplane );
00115     d.snap_to( (void*)( (size_t)HANDLE + 1 ), point );
00116     CPPUNIT_ASSERT_DOUBLES_EQUAL( pplane, point.z(), EPSILON );
00117 }
00118 
00119 void BoundedCylinderDomainTest::test_x_snap_to()
00120 {
00121     const Mesh::EntityHandle HANDLE = 0;
00122     const double XPLANE             = -3;
00123 
00124     // Define a circle in the X=-3 plane, centered at (-3,1,1)
00125     // with a radius of 1
00126     std::vector< Mesh::EntityHandle > handles( 1 );
00127     handles[0] = HANDLE;
00128     BoundedCylinderDomain d( 1, Vector3D( 1, 0, 0 ), Vector3D( 0, 1, 1 ) );
00129     d.create_curve( XPLANE, handles );
00130 
00131     Vector3D point, center( -3, 1, 1 );
00132 
00133     // try a point inside and below the circle
00134     point.set( XPLANE - 2, 1.5, 1.5 );
00135     d.snap_to( HANDLE, point );
00136     CPPUNIT_ASSERT_DOUBLES_EQUAL( XPLANE, point.x(), EPSILON );
00137     CPPUNIT_ASSERT_DOUBLES_EQUAL( point.y(), point.z(), EPSILON );
00138     point -= center;
00139     CPPUNIT_ASSERT_DOUBLES_EQUAL( d.radius(), point.length(), EPSILON );
00140 
00141     // try a point outside and above the circle
00142     point.set( XPLANE + 20, 50, 50 );
00143     d.snap_to( HANDLE, point );
00144     CPPUNIT_ASSERT_DOUBLES_EQUAL( XPLANE, point.x(), EPSILON );
00145     CPPUNIT_ASSERT_DOUBLES_EQUAL( point.y(), point.z(), EPSILON );
00146     point -= center;
00147     CPPUNIT_ASSERT_DOUBLES_EQUAL( d.radius(), point.length(), EPSILON );
00148 
00149     // try a point on the circle
00150     Vector3D on_circ( XPLANE, d.radius(), 0 );
00151     point = on_circ;
00152     d.snap_to( HANDLE, point );
00153     ASSERT_VECTORS_EQUAL( on_circ, point );
00154 
00155     // try a point at the center of the circle
00156     point = center;
00157     d.snap_to( HANDLE, point );
00158     CPPUNIT_ASSERT_DOUBLES_EQUAL( XPLANE, point.x(), EPSILON );
00159     point -= center;
00160     CPPUNIT_ASSERT_DOUBLES_EQUAL( d.radius(), point.length(), EPSILON );
00161 
00162     // try a different handle to make sure a handle not
00163     // on the the circle results in a point not on the circle
00164     double pplane = XPLANE + 5;
00165     point.set( pplane, 0, 0 );
00166     d.snap_to( (void*)( (size_t)HANDLE + 1 ), point );
00167     CPPUNIT_ASSERT_DOUBLES_EQUAL( pplane, point.x(), EPSILON );
00168 }
00169 
00170 // define a VDK file containing four quads quartering a
00171 // square with corners at ( +/-1, +/-1, 0 )
00172 const char vtk_file[] = "# vtk DataFile Version 2.0\n"
00173                         "Mesquite Mesh\n"
00174                         "ASCII\n"
00175                         "DATASET UNSTRUCTURED_GRID\n"
00176                         "POINTS 9 float\n"
00177                         " 0  0 0\n"
00178                         " 1  0 0\n"
00179                         " 1  1 0\n"
00180                         " 0  1 0\n"
00181                         "-1  1 0\n"
00182                         "-1  0 0\n"
00183                         "-1 -1 0\n"
00184                         " 0 -1 0\n"
00185                         " 1 -1 0\n"
00186                         "CELLS 4 20\n"
00187                         "4 0 1 2 3\n"
00188                         "4 0 3 4 5\n"
00189                         "4 0 5 6 7\n"
00190                         "4 0 7 8 1\n"
00191                         "CELL_TYPES 4\n"
00192                         "9 9 9 9\n";
00193 
00194 void BoundedCylinderDomainTest::test_create_curve_from_mesh()
00195 {
00196     // create the input file to read
00197     char filename[] = "BCDTestTMP";
00198     FILE* file      = fopen( filename, "w" );
00199     size_t s        = fwrite( vtk_file, sizeof( vtk_file ), 1, file );
00200     if( !s )
00201     {
00202         remove( filename );
00203         CPPUNIT_ASSERT( s );
00204     }
00205     fclose( file );
00206 
00207     // read the file into a Mesh object
00208     MsqPrintError err( std::cout );
00209     MeshImpl mesh;
00210     mesh.read_vtk( filename, err );
00211     remove( filename );
00212     CPPUNIT_ASSERT( !err );
00213 
00214     // create a domain such that the four corner vertices
00215     // of the 4-quad square mesh lie on the bounding curves
00216     // of the cylinder.
00217     BoundedCylinderDomain d( 1, Vector3D( 0, 1, 0 ) );
00218     d.create_curve( 1, &mesh, EPSILON, err );
00219     CPPUNIT_ASSERT( !err );
00220     d.create_curve( -1, &mesh, EPSILON, err );
00221     CPPUNIT_ASSERT( !err );
00222 
00223     // Get list of vertex coordinates from mesh
00224     std::vector< Mesh::VertexHandle > handles;
00225     mesh.get_all_vertices( handles, err );
00226     CPPUNIT_ASSERT( !err );
00227     CPPUNIT_ASSERT_EQUAL( handles.size(), (size_t)9 );
00228     std::vector< MsqVertex > vertices( handles.size() );
00229     mesh.vertices_get_coordinates( arrptr( handles ), arrptr( vertices ), handles.size(), err );
00230     CPPUNIT_ASSERT( !err );
00231 
00232     // Get the dimension of the geometry each vertex is bound to
00233     std::vector< unsigned short > dims( handles.size() );
00234     d.domain_DoF( arrptr( handles ), arrptr( dims ), handles.size(), err );
00235     CPPUNIT_ASSERT( !err );
00236 
00237     // Check that the four corner vertices are on the curves
00238     // and all others are not
00239     for( size_t i = 0; i < handles.size(); ++i )
00240     {
00241         if( fabs( vertices[i].x() ) < EPSILON || fabs( vertices[i].y() ) < EPSILON )
00242             CPPUNIT_ASSERT_EQUAL( (unsigned short)2, dims[i] );
00243         else
00244             CPPUNIT_ASSERT_EQUAL( (unsigned short)1, dims[i] );
00245     }
00246 }
00247 
00248 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( BoundedCylinderDomainTest, "BoundedCylinderDomainTest" );
00249 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( BoundedCylinderDomainTest, "Unit" );
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines