MOAB: Mesh Oriented datABase  (version 5.2.1)
MsqMeshEntityTest.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: 3; c-tab-always-indent: t; indent-tabs-mode: nil; c-basic-offset: 3
00028 // -*-
00029 //
00030 //   SUMMARY:
00031 //     USAGE:
00032 //
00033 //    AUTHOR: Thomas Leurent <tleurent@mcs.anl.gov>
00034 //       ORG: Argonne National Laboratory
00035 //    E-MAIL: tleurent@mcs.anl.gov
00036 //
00037 // ORIG-DATE: 12-Nov-02 at 18:05:56
00038 //  LAST-MOD:  7-May-03 at 13:28:47 by Thomas Leurent
00039 //
00040 // DESCRIPTION:
00041 // ============
00042 /*! \file MsqMeshEntityTest.cpp
00043 
00044 Unit testing of various functions in the MsqMeshEntity class.
00045 
00046 \author Michael Brewer
00047 \author Thomas Leurent
00048 
00049  */
00050 // DESCRIP-END.
00051 //
00052 
00053 #include "MsqMeshEntity.hpp"
00054 #include "Vector3D.hpp"
00055 #include "PatchData.hpp"
00056 #include "PatchDataInstances.hpp"
00057 #include <math.h>
00058 #include <iostream>
00059 #include <sstream>
00060 #include "UnitUtil.hpp"
00061 
00062 using namespace MBMesquite;
00063 using std::cout;
00064 using std::endl;
00065 
00066 class MsqMeshEntityTest : public CppUnit::TestFixture
00067 {
00068 
00069   private:
00070     CPPUNIT_TEST_SUITE( MsqMeshEntityTest );
00071     CPPUNIT_TEST( test_hex_vertices );
00072     CPPUNIT_TEST( test_centroid_tri );
00073     CPPUNIT_TEST( test_centroid_quad );
00074     CPPUNIT_TEST( test_centroid_hex );
00075     CPPUNIT_TEST( test_unsigned_area );
00076     CPPUNIT_TEST( test_unsigned_area_poly );
00077     CPPUNIT_TEST( test_unsigned_area_tet );
00078     CPPUNIT_TEST( test_unsigned_area_pyr );
00079     CPPUNIT_TEST( test_unsigned_area_pri );
00080     CPPUNIT_TEST( test_unsigned_area_hex );
00081     CPPUNIT_TEST( test_all_nodes );
00082     CPPUNIT_TEST( test_check_element_orientation_linear );
00083     CPPUNIT_TEST( test_check_element_orientation_quadratic );
00084     CPPUNIT_TEST_SUITE_END();
00085 
00086     void test_all_nodes( EntityTopology type, unsigned num_nodes );
00087 
00088   private:
00089     PatchData oneHexPatch;
00090     PatchData oneTetPatch;
00091     PatchData oneQuadPatch;
00092     PatchData oneTriPatch;
00093     Vector3D e1, e2, e3;
00094     double tolEps;
00095 
00096   public:
00097     void setUp()
00098     {
00099         tolEps = 1.e-12;
00100 
00101         // sets up the unit vectors
00102         e1.set( 1, 0, 0 );
00103         e2.set( 0, 1, 0 );
00104         e3.set( 0, 0, 1 );
00105 
00106         MsqPrintError err( cout );
00107 
00108         // creates empty Patch
00109         create_one_hex_patch( oneHexPatch, err );
00110         CPPUNIT_ASSERT( !err );
00111         create_one_tet_patch( oneTetPatch, err );
00112         CPPUNIT_ASSERT( !err );
00113         create_one_tri_patch( oneTriPatch, err );
00114         CPPUNIT_ASSERT( !err );
00115         create_one_quad_patch( oneQuadPatch, err );
00116         CPPUNIT_ASSERT( !err );
00117     }
00118 
00119     void tearDown()
00120     {
00121         destroy_patch_with_domain( oneTriPatch );
00122         destroy_patch_with_domain( oneQuadPatch );
00123     }
00124 
00125   public:
00126     MsqMeshEntityTest() {}
00127 
00128     void test_hex_vertices()
00129     {
00130         MsqPrintError err( cout );
00131         // prints out the vertices.
00132         const MsqVertex* ideal_vertices = oneHexPatch.get_vertex_array( err );
00133         CPPUNIT_ASSERT( !err );
00134         size_t num_vtx = oneHexPatch.num_nodes();
00135         CPPUNIT_ASSERT_EQUAL( size_t( 8 ), num_vtx );
00136 
00137         MsqVertex vtx;
00138 
00139         vtx.set( 1, 1, 1 );
00140         CPPUNIT_ASSERT_EQUAL( vtx, ideal_vertices[0] );
00141 
00142         vtx.set( 2, 2, 2 );
00143         CPPUNIT_ASSERT_EQUAL( vtx, ideal_vertices[6] );
00144 
00145         vtx.set( 1, 2, 2 );
00146         CPPUNIT_ASSERT_EQUAL( vtx, ideal_vertices[7] );
00147     }
00148 
00149     //! test the centroid of the first element in the Patch
00150     void test_centroid( PatchData& pd, Vector3D& correct )
00151     {
00152         MsqPrintError err( cout );
00153         double eps = 1e-6;
00154         Vector3D centroid;
00155 
00156         MsqMeshEntity* elem = pd.get_element_array( err );
00157         CPPUNIT_ASSERT( !err );
00158         elem->get_centroid( centroid, pd, err );
00159         CPPUNIT_ASSERT( !err );
00160 
00161         //     cout << "centroid: "<< centroid <<endl;
00162         //     cout << "correct: "<< correct <<endl;
00163 
00164         for( int i = 0; i < 3; ++i )
00165             CPPUNIT_ASSERT_DOUBLES_EQUAL( centroid[i], correct[i], eps );
00166     }
00167 
00168     void test_centroid_tri()
00169     {
00170         Vector3D correct( 1.5, 1 + 1 / ( 2.0 * sqrt( 3.0 ) ), 1.0 );
00171         test_centroid( oneTriPatch, correct );
00172     }
00173 
00174     void test_centroid_quad()
00175     {
00176         Vector3D correct( 1.5, 1.5, 1.0 );
00177         test_centroid( oneQuadPatch, correct );
00178     }
00179 
00180     void test_centroid_hex()
00181     {
00182         Vector3D correct( 1.5, 1.5, 1.5 );
00183         test_centroid( oneHexPatch, correct );
00184     }
00185 
00186     void test_unsigned_area()
00187     {
00188         MsqPrintError err( cout );
00189         MsqMeshEntity* tri = oneTriPatch.get_element_array( err );
00190         CPPUNIT_ASSERT( !err );
00191         CPPUNIT_ASSERT( fabs( tri->compute_unsigned_area( oneTriPatch, err ) - ( sqrt( 3.0 ) / 4.0 ) ) < tolEps );
00192         MsqMeshEntity* quad = oneQuadPatch.get_element_array( err );
00193         CPPUNIT_ASSERT( !err );
00194         CPPUNIT_ASSERT( fabs( quad->compute_unsigned_area( oneQuadPatch, err ) - 1.0 ) < tolEps );
00195     }
00196 
00197     void test_unsigned_area_poly();
00198     void test_unsigned_area_tet();
00199     void test_unsigned_area_pyr();
00200     void test_unsigned_area_pri();
00201     void test_unsigned_area_hex();
00202     void test_all_nodes();
00203 
00204     void test_unsigned_area_common( EntityTopology type, const double* coords, double expected );
00205 
00206     void test_check_element_orientation_linear();
00207     void test_check_element_orientation_quadratic();
00208     void test_check_element_orientation( EntityTopology type, int nodes );
00209 };
00210 
00211 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MsqMeshEntityTest, "MsqMeshEntityTest" );
00212 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MsqMeshEntityTest, "Unit" );
00213 
00214 const size_t conn[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
00215 const bool fixed[]  = { false, false, false, false, false, false, false, false };
00216 
00217 void MsqMeshEntityTest::test_unsigned_area_poly()
00218 {
00219     const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0.5, 1.5, 0, 0, 1, 0 };
00220     size_t n_vtx          = 5;
00221     EntityTopology type   = POLYGON;
00222     MsqError err;
00223 
00224     PatchData pd;
00225     pd.fill( n_vtx, coords, 1, &type, &n_vtx, conn, fixed, err );
00226     ASSERT_NO_ERROR( err );
00227 
00228     double a = pd.element_by_index( 0 ).compute_unsigned_area( pd, err );
00229     ASSERT_NO_ERROR( err );
00230     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.25, a, 1e-8 );
00231 }
00232 
00233 void MsqMeshEntityTest::test_unsigned_area_common( EntityTopology type, const double* coords, double expected )
00234 {
00235     MsqError err;
00236     PatchData pd;
00237 
00238     pd.fill( TopologyInfo::corners( type ), coords, 1, type, conn, fixed, err );
00239     ASSERT_NO_ERROR( err );
00240 
00241     double a = pd.element_by_index( 0 ).compute_unsigned_area( pd, err );
00242     ASSERT_NO_ERROR( err );
00243 
00244     CPPUNIT_ASSERT_DOUBLES_EQUAL( expected, a, 1e-8 );
00245 }
00246 
00247 void MsqMeshEntityTest::test_unsigned_area_tet()
00248 {
00249     const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 };
00250     test_unsigned_area_common( TETRAHEDRON, coords, 1.0 / 6.0 );
00251 }
00252 
00253 void MsqMeshEntityTest::test_unsigned_area_pyr()
00254 {
00255     const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1 };
00256     test_unsigned_area_common( PYRAMID, coords, 1.0 / 3.0 );
00257 
00258     const double pyr_coords[] = { -1, -1, -1, 1, -1, -1, 1, 1, -1, -1, 1, -1, 0, 0, 0 };
00259     test_unsigned_area_common( PYRAMID, pyr_coords, 4.0 / 3.0 );
00260 }
00261 
00262 void MsqMeshEntityTest::test_unsigned_area_pri()
00263 {
00264     const double coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 0, 1, 1 };
00265     test_unsigned_area_common( PRISM, coords, 0.5 );
00266 
00267     const double tet_coords[] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1 };
00268     test_unsigned_area_common( PRISM, tet_coords, 1.0 / 6.0 );
00269 }
00270 
00271 void MsqMeshEntityTest::test_unsigned_area_hex()
00272 {
00273     const double coords[] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1 };
00274     test_unsigned_area_common( HEXAHEDRON, coords, 1.0 );
00275 
00276     const double coords2[] = { 0, 0, 0, 2, 0, 0, 2, 2, 0, 0, 2, 0, 0, 0, 2, 2, 0, 2, 2, 2, 2, 0, 2, 2 };
00277     test_unsigned_area_common( HEXAHEDRON, coords2, 8.0 );
00278 
00279     const double pyr_coords[] = { -1, -1, 0, 1, -1, 0, 1, 1, 0, -1, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 1 };
00280     test_unsigned_area_common( HEXAHEDRON, pyr_coords, 4.0 / 3.0 );
00281 }
00282 
00283 void MsqMeshEntityTest::test_all_nodes( EntityTopology type, unsigned num_nodes )
00284 {
00285     const unsigned num_vtx     = 27;
00286     double coords[3 * num_vtx] = { 0.0 };
00287     size_t conn[num_vtx];
00288     for( size_t i = 0; i < num_vtx; ++i )
00289         conn[i] = i;
00290     bool fixed[num_vtx] = { false };
00291     CPPUNIT_ASSERT( num_nodes <= num_vtx );
00292 
00293     MsqError err;
00294     PatchData pd;
00295     size_t n = num_nodes;
00296     pd.fill( num_nodes, coords, 1, &type, &n, conn, fixed, err );
00297     ASSERT_NO_ERROR( err );
00298 
00299     MsqMeshEntity& elem = pd.element_by_index( 0 );
00300     NodeSet all         = elem.all_nodes( err );
00301     ASSERT_NO_ERROR( err );
00302     CPPUNIT_ASSERT_EQUAL( num_nodes, all.num_nodes() );
00303     CPPUNIT_ASSERT( all.have_any_corner_node() );
00304     bool mid_edge, mid_face, mid_reg;
00305     TopologyInfo::higher_order( type, num_nodes, mid_edge, mid_face, mid_reg, err );
00306     ASSERT_NO_ERROR( err );
00307     CPPUNIT_ASSERT_EQUAL( mid_edge, !!all.have_any_mid_edge_node() );
00308     CPPUNIT_ASSERT_EQUAL( mid_face, !!all.have_any_mid_face_node() );
00309     CPPUNIT_ASSERT_EQUAL( mid_reg, !!all.have_any_mid_region_node() );
00310 }
00311 
00312 void MsqMeshEntityTest::test_all_nodes()
00313 {
00314     test_all_nodes( TRIANGLE, 3 );
00315     test_all_nodes( TRIANGLE, 4 );
00316     test_all_nodes( TRIANGLE, 6 );
00317     test_all_nodes( TRIANGLE, 7 );
00318 
00319     test_all_nodes( QUADRILATERAL, 4 );
00320     test_all_nodes( QUADRILATERAL, 5 );
00321     test_all_nodes( QUADRILATERAL, 8 );
00322     test_all_nodes( QUADRILATERAL, 9 );
00323 
00324     test_all_nodes( TETRAHEDRON, 4 );
00325     test_all_nodes( TETRAHEDRON, 5 );
00326     test_all_nodes( TETRAHEDRON, 10 );
00327     test_all_nodes( TETRAHEDRON, 11 );
00328     test_all_nodes( TETRAHEDRON, 8 );
00329     test_all_nodes( TETRAHEDRON, 9 );
00330     test_all_nodes( TETRAHEDRON, 14 );
00331     test_all_nodes( TETRAHEDRON, 15 );
00332 
00333     test_all_nodes( HEXAHEDRON, 8 );
00334     test_all_nodes( HEXAHEDRON, 9 );
00335     test_all_nodes( HEXAHEDRON, 20 );
00336     test_all_nodes( HEXAHEDRON, 21 );
00337     test_all_nodes( HEXAHEDRON, 14 );
00338     test_all_nodes( HEXAHEDRON, 15 );
00339     test_all_nodes( HEXAHEDRON, 26 );
00340     test_all_nodes( HEXAHEDRON, 27 );
00341 }
00342 
00343 void MsqMeshEntityTest::test_check_element_orientation_linear()
00344 {
00345     const EntityTopology types[] = { TRIANGLE, QUADRILATERAL, TETRAHEDRON, PYRAMID, PRISM, HEXAHEDRON };
00346     const int num_types          = sizeof( types ) / sizeof( types[0] );
00347 
00348     for( int i = 0; i < num_types; ++i )
00349     {
00350         test_check_element_orientation( types[i], TopologyInfo::corners( types[i] ) );
00351     }
00352 }
00353 
00354 void MsqMeshEntityTest::test_check_element_orientation_quadratic()
00355 {
00356     struct ElemType
00357     {
00358         EntityTopology topo;
00359         unsigned nodes;
00360     };
00361 
00362     const ElemType types[] = { { TRIANGLE, 6 }, { QUADRILATERAL, 8 }, { QUADRILATERAL, 9 }, { TETRAHEDRON, 10 } };
00363     const int num_types    = sizeof( types ) / sizeof( types[0] );
00364 
00365     for( int i = 0; i < num_types; ++i )
00366     {
00367         test_check_element_orientation( types[i].topo, types[i].nodes );
00368     }
00369 }
00370 
00371 void MsqMeshEntityTest::test_check_element_orientation( EntityTopology type, int nodes )
00372 {
00373     // get an ideal element
00374     MsqError err;
00375     PatchData pd;
00376     create_ideal_element_patch( pd, type, nodes, err );
00377     ASSERT_NO_ERROR( err );
00378     CPPUNIT_ASSERT_EQUAL( (size_t)1, pd.num_elements() );
00379     CPPUNIT_ASSERT_EQUAL( (size_t)nodes, pd.num_nodes() );
00380     MsqMeshEntity& elem = pd.element_by_index( 0 );
00381     CPPUNIT_ASSERT_EQUAL( (size_t)nodes, elem.node_count() );
00382     CPPUNIT_ASSERT_EQUAL( type, elem.get_element_type() );
00383     const size_t* conn = elem.get_vertex_index_array();
00384 
00385     // test that ideal element is not reported as inverted
00386     int inverted, tested;
00387     elem.check_element_orientation( pd, inverted, tested, err );
00388     ASSERT_NO_ERROR( err );
00389     CPPUNIT_ASSERT_EQUAL( 0, inverted );
00390     CPPUNIT_ASSERT( tested > 0 );
00391 
00392     bool mids[4] = { false };
00393     TopologyInfo::higher_order( type, nodes, mids[1], mids[2], mids[3], err );MSQ_ERRRTN( err );
00394 
00395     // invert element at each vertex and test
00396     Vector3D centroid;
00397     elem.get_centroid( centroid, pd, err );
00398     ASSERT_NO_ERROR( err );
00399     for( int i = 0; i < nodes; ++i )
00400     {
00401         unsigned dim, num;
00402         TopologyInfo::side_from_higher_order( type, nodes, i, dim, num, err );
00403         ASSERT_NO_ERROR( err );
00404         const Vector3D old_pos = pd.vertex_by_index( conn[i] );
00405         Vector3D new_pos       = old_pos;
00406         if( dim == TopologyInfo::dimension( type ) )
00407         {
00408             // move mid-element node 3/4 of the way to corner 0
00409             new_pos += 3 * pd.vertex_by_index( conn[0] );
00410             new_pos *= 0.25;
00411         }
00412         else if( dim == 0 )
00413         {  // if a corner vertex
00414             if( type == TRIANGLE || type == TETRAHEDRON )
00415             {
00416                 // move tri/tet vertex past opposite side of element
00417                 new_pos += 2 * ( centroid - old_pos );
00418             }
00419             else if( mids[1] )
00420             {
00421                 // if have mid-edge nodes move 3/4 of the way to center vertex
00422                 new_pos += 3 * centroid;
00423                 new_pos *= 0.25;
00424             }
00425             else
00426             {
00427                 // move vertex past centroid
00428                 new_pos += 1.5 * ( centroid - old_pos );
00429             }
00430         }
00431         else
00432         {
00433             // otherwise move vertex past centroid
00434             new_pos += 2.5 * ( centroid - old_pos );
00435         }
00436 
00437         pd.set_vertex_coordinates( new_pos, conn[i], err );
00438         ASSERT_NO_ERROR( err );
00439 
00440         // test that element is inverted
00441         inverted = tested = 0;
00442         elem.check_element_orientation( pd, inverted, tested, err );
00443         ASSERT_NO_ERROR( err );
00444         std::ostringstream str;
00445         str << TopologyInfo::short_name( type ) << nodes << " Vertex " << i << " (Dimension " << dim << " Index " << num
00446             << ")";
00447         CppUnit::Message m( "MsqMeshEntity failed to detect inverted element" );
00448         m.addDetail( str.str() );
00449         ASSERT_MESSAGE( m, inverted > 0 );
00450 
00451         // move vertex back to ideal position
00452         pd.set_vertex_coordinates( old_pos, conn[i], err );
00453         ASSERT_NO_ERROR( err );
00454     }
00455 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines