MOAB: Mesh Oriented datABase  (version 5.3.1)
IdealElementTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Sandia National Laboratories.  Developed at the
00005     University of Wisconsin--Madison under SNL contract number
00006     624796.  The U.S. Government and the University of Wisconsin
00007     retian certain rights to 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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file IdealElementTest.cpp
00028  *  \brief Test the ExtraData functionality of the PatchData class
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "IdealElements.hpp"
00034 #include "TopologyInfo.hpp"
00035 #include "MsqError.hpp"
00036 #include "Vector3D.hpp"
00037 #include "cppunit/extensions/HelperMacros.h"
00038 
00039 using namespace MBMesquite;
00040 
00041 class IdealElementTest : public CppUnit::TestFixture
00042 {
00043   private:
00044     CPPUNIT_TEST_SUITE( IdealElementTest );
00045     CPPUNIT_TEST( test_unit_tri );
00046     CPPUNIT_TEST( test_unit_quad );
00047     CPPUNIT_TEST( test_unit_tet );
00048     CPPUNIT_TEST( test_unit_pyr );
00049     CPPUNIT_TEST( test_unit_wdg );
00050     CPPUNIT_TEST( test_unit_hex );
00051     CPPUNIT_TEST( test_side_height_pyr );
00052     CPPUNIT_TEST( test_unit_edge_tri );
00053     CPPUNIT_TEST( test_unit_edge_quad );
00054     CPPUNIT_TEST( test_unit_edge_tet );
00055     CPPUNIT_TEST( test_unit_edge_pyr );
00056     CPPUNIT_TEST( test_unit_edge_wdg );
00057     CPPUNIT_TEST( test_unit_edge_hex );
00058     CPPUNIT_TEST( test_unit_height_pyr );
00059     CPPUNIT_TEST_SUITE_END();
00060 
00061   public:
00062     void test_unit_tri();
00063     void test_unit_quad();
00064     void test_unit_tet();
00065     void test_unit_pyr();
00066     void test_unit_wdg();
00067     void test_unit_hex();
00068     void test_unit_edge_tri();
00069     void test_unit_edge_quad();
00070     void test_unit_edge_tet();
00071     void test_unit_edge_pyr();
00072     void test_unit_edge_wdg();
00073     void test_unit_edge_hex();
00074     void test_side_height_pyr();
00075     void test_unit_height_pyr();
00076 };
00077 
00078 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( IdealElementTest, "IdealElementTest" );
00079 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( IdealElementTest, "Unit" );
00080 
00081 static double distance_from_origin( unsigned n, const Vector3D* vtx );
00082 static bool unit_edge_lengths( EntityTopology type, const Vector3D* coords );
00083 static bool equal_edge_lengths( EntityTopology type, const Vector3D* coords );
00084 
00085 void IdealElementTest::test_unit_tri()
00086 {
00087     const Vector3D* coords = unit_element( TRIANGLE );
00088     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 3, coords ), 1e-6 );
00089     CPPUNIT_ASSERT( equal_edge_lengths( TRIANGLE, coords ) );
00090 
00091     Vector3D v1 = coords[1] - coords[0];
00092     Vector3D v2 = coords[2] - coords[0];
00093     double area = 0.5 * ( v1 * v2 ).length();
00094     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, area, 1e-6 );
00095 }
00096 
00097 void IdealElementTest::test_unit_quad()
00098 {
00099     const Vector3D* coords = unit_element( QUADRILATERAL );
00100     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 4, coords ), 1e-6 );
00101     CPPUNIT_ASSERT( equal_edge_lengths( QUADRILATERAL, coords ) );
00102 
00103     Vector3D p1 = 0.5 * ( coords[0] + coords[1] );
00104     Vector3D p2 = 0.5 * ( coords[1] + coords[2] );
00105     Vector3D p3 = 0.5 * ( coords[2] + coords[3] );
00106     Vector3D p4 = 0.5 * ( coords[3] + coords[0] );
00107     Vector3D v1 = p1 - p3;
00108     Vector3D v2 = p2 - p4;
00109     double area = ( v1 * v2 ).length();
00110     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, area, 1e-6 );
00111 }
00112 
00113 void IdealElementTest::test_unit_tet()
00114 {
00115     const Vector3D* coords = unit_element( TETRAHEDRON );
00116     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 4, coords ), 1e-6 );
00117     CPPUNIT_ASSERT( equal_edge_lengths( TETRAHEDRON, coords ) );
00118 
00119     Vector3D v1 = coords[1] - coords[0];
00120     Vector3D v2 = coords[2] - coords[0];
00121     Vector3D v3 = coords[3] - coords[0];
00122     double vol  = ( v3 % ( v1 * v2 ) ) / 6.0;
00123     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, vol, 1e-6 );
00124 }
00125 
00126 void IdealElementTest::test_unit_pyr()
00127 {
00128     const Vector3D* coords = unit_element( PYRAMID );
00129     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 5, coords ), 1e-6 );
00130     CPPUNIT_ASSERT( equal_edge_lengths( PYRAMID, coords ) );
00131 
00132     Vector3D p1 = 0.5 * ( coords[0] + coords[1] );
00133     Vector3D p2 = 0.5 * ( coords[1] + coords[2] );
00134     Vector3D p3 = 0.5 * ( coords[2] + coords[3] );
00135     Vector3D p4 = 0.5 * ( coords[3] + coords[0] );
00136     Vector3D v1 = p1 - p3;
00137     Vector3D v2 = p2 - p4;
00138     double area = ( v1 * v2 ).length();
00139 
00140     Vector3D n = ( v1 * v2 );
00141     n.normalize();
00142     Vector3D c    = 0.25 * ( coords[0] + coords[1] + coords[2] + coords[3] );
00143     double height = n % ( coords[4] - c );
00144     double vol    = area * height / 3.0;
00145     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, vol, 1e-6 );
00146 }
00147 
00148 void IdealElementTest::test_unit_wdg()
00149 {
00150     const Vector3D* coords = unit_element( PRISM );
00151     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 6, coords ), 1e-6 );
00152     CPPUNIT_ASSERT( equal_edge_lengths( PRISM, coords ) );
00153 
00154     Vector3D v1 = coords[1] - coords[0];
00155     Vector3D v2 = coords[2] - coords[0];
00156     Vector3D v3 = coords[3] - coords[0];
00157     double vol  = 0.5 * v3 % ( v1 * v2 );
00158     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, vol, 1e-6 );
00159 }
00160 
00161 void IdealElementTest::test_unit_hex()
00162 {
00163     const Vector3D* coords = unit_element( HEXAHEDRON );
00164     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 8, coords ), 1e-6 );
00165     CPPUNIT_ASSERT( equal_edge_lengths( HEXAHEDRON, coords ) );
00166 
00167     Vector3D p1 = 0.25 * ( coords[0] + coords[1] + coords[2] + coords[3] );
00168     Vector3D p2 = 0.25 * ( coords[4] + coords[5] + coords[6] + coords[7] );
00169     Vector3D p3 = 0.25 * ( coords[0] + coords[1] + coords[5] + coords[4] );
00170     Vector3D p4 = 0.25 * ( coords[7] + coords[6] + coords[6] + coords[7] );
00171     Vector3D p5 = 0.25 * ( coords[0] + coords[3] + coords[4] + coords[7] );
00172     Vector3D p6 = 0.25 * ( coords[1] + coords[2] + coords[6] + coords[5] );
00173     Vector3D v1 = p1 - p2;
00174     Vector3D v2 = p3 - p4;
00175     Vector3D v3 = p5 - p6;
00176     double vol  = v3 % ( v1 * v2 );
00177     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, vol, 1e-6 );
00178 }
00179 
00180 void IdealElementTest::test_side_height_pyr()
00181 {
00182     const Vector3D* coords = unit_element( PYRAMID, true );
00183     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 5, coords ), 1e-6 );
00184     CPPUNIT_ASSERT( equal_edge_lengths( QUADRILATERAL, coords ) );
00185 
00186     Vector3D p1 = 0.5 * ( coords[0] + coords[1] );
00187     Vector3D p2 = 0.5 * ( coords[1] + coords[2] );
00188     Vector3D p3 = 0.5 * ( coords[2] + coords[3] );
00189     Vector3D p4 = 0.5 * ( coords[3] + coords[0] );
00190     Vector3D v1 = p1 - p3;
00191     Vector3D v2 = p2 - p4;
00192     double area = ( v1 * v2 ).length();
00193 
00194     Vector3D n = ( v1 * v2 );
00195     n.normalize();
00196     Vector3D c    = 0.25 * ( coords[0] + coords[1] + coords[2] + coords[3] );
00197     double height = n % ( coords[4] - c );
00198     double vol    = area * height / 3.0;
00199     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, vol, 1e-6 );
00200     CPPUNIT_ASSERT_DOUBLES_EQUAL( height, ( coords[0] - coords[1] ).length(), 1e-6 );
00201 }
00202 
00203 void IdealElementTest::test_unit_edge_tri()
00204 {
00205     const Vector3D* coords = unit_edge_element( TRIANGLE );
00206     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 3, coords ), 1e-6 );
00207     CPPUNIT_ASSERT( unit_edge_lengths( TRIANGLE, coords ) );
00208 }
00209 
00210 void IdealElementTest::test_unit_edge_quad()
00211 {
00212     const Vector3D* coords = unit_edge_element( QUADRILATERAL );
00213     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 4, coords ), 1e-6 );
00214     CPPUNIT_ASSERT( unit_edge_lengths( QUADRILATERAL, coords ) );
00215 }
00216 
00217 void IdealElementTest::test_unit_edge_tet()
00218 {
00219     const Vector3D* coords = unit_edge_element( TETRAHEDRON );
00220     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 4, coords ), 1e-6 );
00221     CPPUNIT_ASSERT( unit_edge_lengths( TETRAHEDRON, coords ) );
00222 }
00223 
00224 void IdealElementTest::test_unit_edge_pyr()
00225 {
00226     const Vector3D* coords = unit_edge_element( PYRAMID );
00227     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 5, coords ), 1e-6 );
00228     CPPUNIT_ASSERT( unit_edge_lengths( PYRAMID, coords ) );
00229 }
00230 
00231 void IdealElementTest::test_unit_edge_wdg()
00232 {
00233     const Vector3D* coords = unit_edge_element( PRISM );
00234     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 6, coords ), 1e-6 );
00235     CPPUNIT_ASSERT( unit_edge_lengths( PRISM, coords ) );
00236 }
00237 
00238 void IdealElementTest::test_unit_edge_hex()
00239 {
00240     const Vector3D* coords = unit_edge_element( HEXAHEDRON );
00241     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 8, coords ), 1e-6 );
00242     CPPUNIT_ASSERT( unit_edge_lengths( HEXAHEDRON, coords ) );
00243 }
00244 
00245 void IdealElementTest::test_unit_height_pyr()
00246 {
00247     const Vector3D* coords = unit_edge_element( PYRAMID, true );
00248     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, distance_from_origin( 5, coords ), 1e-6 );
00249     CPPUNIT_ASSERT( unit_edge_lengths( QUADRILATERAL, coords ) );
00250 
00251     Vector3D p1 = 0.5 * ( coords[0] + coords[1] );
00252     Vector3D p2 = 0.5 * ( coords[1] + coords[2] );
00253     Vector3D p3 = 0.5 * ( coords[2] + coords[3] );
00254     Vector3D p4 = 0.5 * ( coords[3] + coords[0] );
00255     Vector3D v1 = p1 - p3;
00256     Vector3D v2 = p2 - p4;
00257 
00258     Vector3D n = ( v1 * v2 );
00259     n.normalize();
00260     Vector3D c    = 0.25 * ( coords[0] + coords[1] + coords[2] + coords[3] );
00261     double height = n % ( coords[4] - c );
00262     CPPUNIT_ASSERT_DOUBLES_EQUAL( height, ( coords[0] - coords[1] ).length(), 1e-6 );
00263 }
00264 
00265 static void get_edge_lengths( EntityTopology type, const Vector3D* coords, double& min, double& max )
00266 {
00267     MsqError err;
00268     unsigned n = TopologyInfo::edges( type );
00269     min        = HUGE_VAL;
00270     max        = -HUGE_VAL;
00271     for( unsigned j = 0; j < n; ++j )
00272     {
00273         const unsigned* edge = TopologyInfo::edge_vertices( type, j, err );
00274         CPPUNIT_ASSERT( !err );
00275         double len = ( coords[edge[0]] - coords[edge[1]] ).length();
00276         if( min > len ) min = len;
00277         if( max < len ) max = len;
00278     }
00279 }
00280 
00281 static bool unit_edge_lengths( EntityTopology type, const Vector3D* coords )
00282 {
00283     double min, max;
00284     get_edge_lengths( type, coords, min, max );
00285     return fabs( 1.0 - min ) < 1e-6 && fabs( 1.0 - max ) < 1e-6;
00286 }
00287 
00288 static bool equal_edge_lengths( EntityTopology type, const Vector3D* coords )
00289 {
00290     double min, max;
00291     get_edge_lengths( type, coords, min, max );
00292     return ( max - min ) < 1e-6;
00293 }
00294 
00295 static double distance_from_origin( unsigned n, const Vector3D* vtx )
00296 {
00297     Vector3D sum( 0, 0, 0 );
00298     for( unsigned i = 0; i < n; ++i )
00299         sum += vtx[i];
00300     if( n == 5 )  // pyrmid centroid is not mean of corners
00301         sum += ( vtx[4] / 3.0 );
00302 
00303     return sum.length();
00304 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines