MOAB: Mesh Oriented datABase  (version 5.4.1)
MappingFunctionTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2009 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     retain 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     (2009) [email protected]
00024 
00025   ***************************************************************** */
00026 
00027 /** \file MappingFunctionTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "UnitUtil.hpp"
00034 #include "LinearTetrahedron.hpp"
00035 #include "LinearTriangle.hpp"
00036 #include "PatchData.hpp"
00037 #include "IdealElements.hpp"
00038 #include <algorithm>
00039 
00040 using namespace MBMesquite;
00041 
00042 class MappingFunctionTest : public CppUnit::TestFixture
00043 {
00044   private:
00045     CPPUNIT_TEST_SUITE( MappingFunctionTest );
00046     CPPUNIT_TEST( test_jacobian_2d );
00047     CPPUNIT_TEST( test_jacobian_3d );
00048     CPPUNIT_TEST( test_ideal_2d );
00049     CPPUNIT_TEST( test_ideal_3d );
00050     CPPUNIT_TEST_SUITE_END();
00051 
00052     LinearTriangle tri;
00053     LinearTetrahedron tet;
00054 
00055   public:
00056     void test_jacobian_2d();
00057     void test_jacobian_3d();
00058     void test_ideal_2d();
00059     void test_ideal_3d();
00060 
00061     MsqMatrix< 3, 2 > calculate_jacobian_2d( Sample s, const Vector3D* coords );
00062     MsqMatrix< 3, 3 > calculate_jacobian_3d( Sample s, const Vector3D* coords );
00063 };
00064 
00065 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MappingFunctionTest, "MappingFunctionTest" );
00066 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( MappingFunctionTest, "Unit" );
00067 
00068 void MappingFunctionTest::test_jacobian_2d()
00069 {
00070     MsqError err;
00071     PatchData pd;
00072     const size_t conn[3]      = { 0, 1, 2 };
00073     const bool fixed[3]       = { false, false, false };
00074     const double coords[3][3] = { { 0.5, 0.5, 0.0 }, { 2.0, 0.2, 0.0 }, { 0.7, 0.8, 0.0 } };
00075 
00076     size_t indices_exp[3], indices_act[3], num_exp, num_act;
00077     MsqVector< 2 > coeff_exp[3], coeff_act[3];
00078     MsqMatrix< 3, 2 > J_exp, J_act;
00079     for( int r = 0; r < 3; ++r )
00080     {
00081         const double verts[] = { coords[r][0],
00082                                  coords[r][1],
00083                                  coords[r][2],
00084                                  coords[( r + 1 ) % 3][0],
00085                                  coords[( r + 1 ) % 3][1],
00086                                  coords[( r + 1 ) % 3][2],
00087                                  coords[( r + 2 ) % 3][0],
00088                                  coords[( r + 2 ) % 3][1],
00089                                  coords[( r + 2 ) % 3][2] };
00090         pd.fill( 3, verts, 1, TRIANGLE, conn, fixed, err );
00091         ASSERT_NO_ERROR( err );
00092 
00093         for( int d = 0; d < 3; ++d )
00094         {
00095             int n = ( d == 2 ) ? 1 : 3;
00096             for( int s = 0; s < n; ++s )
00097             {
00098                 tri.MappingFunction2D::jacobian( pd, 0, NodeSet(), Sample( d, s ), indices_act, coeff_act, num_act,
00099                                                  J_act, err );
00100                 ASSERT_NO_ERROR( err );
00101                 CPPUNIT_ASSERT( num_act <= 3 );
00102                 CPPUNIT_ASSERT( num_act > 0 );
00103 
00104                 tri.derivatives( Sample( d, s ), NodeSet(), indices_exp, coeff_exp, num_exp, err );
00105                 ASSERT_NO_ERROR( err );
00106                 CPPUNIT_ASSERT_EQUAL( num_exp, num_act );
00107 
00108                 J_exp = MsqMatrix< 3, 2 >( 0.0 );  // zero
00109                 for( size_t j = 0; j < num_exp; ++j )
00110                 {
00111                     size_t idx = pd.element_by_index( 0 ).get_vertex_index_array()[j];
00112                     size_t k   = std::find( indices_act, indices_act + num_act, idx ) - indices_act;
00113                     CPPUNIT_ASSERT( k < num_act );  // found
00114                     ASSERT_MATRICES_EQUAL( coeff_exp[j], coeff_act[k], 1e-10 );
00115 
00116                     J_exp += MsqVector< 3 >( pd.vertex_by_index( idx ).to_array() ) * transpose( coeff_exp[j] );
00117                 }
00118 
00119                 ASSERT_MATRICES_EQUAL( J_exp, J_act, 1e-6 );
00120             }
00121         }
00122     }
00123 }
00124 
00125 void MappingFunctionTest::test_jacobian_3d()
00126 {
00127     MsqError err;
00128     PatchData pd;
00129     const size_t conn[4]      = { 0, 1, 2, 3 };
00130     const bool fixed[4]       = { false, false, false, false };
00131     const double coords[4][3] = { { 0.5, 0.5, 0.0 }, { 2.0, 0.2, 0.0 }, { 0.7, 0.8, 0.0 }, { 0.6, 1.1, 13. } };
00132 
00133     size_t indices_exp[4], indices_act[4], num_exp, num_act;
00134     MsqVector< 3 > coeff_exp[4], coeff_act[4];
00135     MsqMatrix< 3, 3 > J_exp, J_act;
00136     for( int r = 0; r < 3; ++r )
00137     {
00138         const double verts[] = { coords[r][0],
00139                                  coords[r][1],
00140                                  coords[r][2],
00141                                  coords[( r + 1 ) % 3][0],
00142                                  coords[( r + 1 ) % 3][1],
00143                                  coords[( r + 1 ) % 3][2],
00144                                  coords[( r + 2 ) % 3][0],
00145                                  coords[( r + 2 ) % 3][1],
00146                                  coords[( r + 2 ) % 3][2],
00147                                  coords[3][0],
00148                                  coords[3][1],
00149                                  coords[3][2] };
00150         pd.fill( 4, verts, 1, TETRAHEDRON, conn, fixed, err );
00151         ASSERT_NO_ERROR( err );
00152 
00153         for( int d = 0; d < 3; ++d )
00154         {
00155             int n = ( d == 3 ) ? 1 : ( d == 1 ) ? 6 : 4;
00156             for( int s = 0; s < n; ++s )
00157             {
00158                 tet.MappingFunction3D::jacobian( pd, 0, NodeSet(), Sample( d, s ), indices_act, coeff_act, num_act,
00159                                                  J_act, err );
00160                 ASSERT_NO_ERROR( err );
00161                 CPPUNIT_ASSERT( num_act <= 4 );
00162                 CPPUNIT_ASSERT( num_act > 0 );
00163 
00164                 tet.derivatives( Sample( d, s ), NodeSet(), indices_exp, coeff_exp, num_exp, err );
00165                 ASSERT_NO_ERROR( err );
00166                 CPPUNIT_ASSERT_EQUAL( num_exp, num_act );
00167 
00168                 J_exp = MsqMatrix< 3, 3 >( 0.0 );  // zero
00169                 for( size_t j = 0; j < num_exp; ++j )
00170                 {
00171                     size_t idx = pd.element_by_index( 0 ).get_vertex_index_array()[j];
00172                     size_t k   = std::find( indices_act, indices_act + num_act, idx ) - indices_act;
00173                     CPPUNIT_ASSERT( k < num_act );  // found
00174                     ASSERT_MATRICES_EQUAL( coeff_exp[j], coeff_act[k], 1e-10 );
00175 
00176                     J_exp += MsqVector< 3 >( pd.vertex_by_index( idx ).to_array() ) * transpose( coeff_exp[j] );
00177                 }
00178 
00179                 ASSERT_MATRICES_EQUAL( J_exp, J_act, 1e-6 );
00180             }
00181         }
00182     }
00183 }
00184 
00185 void MappingFunctionTest::test_ideal_2d()
00186 {
00187     MsqError err;
00188     MsqMatrix< 3, 2 > W_prime;
00189     tri.MappingFunction2D::ideal( Sample( 2, 0 ), W_prime, err );
00190     ASSERT_NO_ERROR( err );
00191 
00192     // for this test that everything is in the xy-plane
00193     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 0 ), 1e-12 );
00194     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 1 ), 1e-12 );
00195     MsqMatrix< 2, 2 > W( W_prime.data() );
00196     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( W ), 1e-6 );
00197 
00198     const Vector3D* verts = unit_edge_element( TRIANGLE );
00199     CPPUNIT_ASSERT( verts );
00200 
00201     size_t indices[3], num;
00202     MsqVector< 2 > coeff[3];
00203     tri.derivatives( Sample( 2, 0 ), NodeSet(), indices, coeff, num, err );
00204     ASSERT_NO_ERROR( err );
00205 
00206     MsqMatrix< 3, 2 > J_exp( 0.0 );  // zero
00207     for( size_t j = 0; j < num; ++j )
00208         J_exp += MsqVector< 3 >( verts[j].to_array() ) * transpose( coeff[j] );
00209 
00210     // for this test that everything is in the xy-plane
00211     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 0 ), 1e-12 );
00212     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 1 ), 1e-12 );
00213     MsqMatrix< 2, 2 > W_exp( J_exp.data() );
00214     W_exp /= sqrt( det( W_exp ) );
00215 
00216     // Matrices should be a rotation of each other.
00217     // First, calculate tentative rotation matrix
00218     MsqMatrix< 2, 2 > R = inverse( W_exp ) * W;
00219     // next check that it is a rotation
00220     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00221     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00222 }
00223 
00224 void MappingFunctionTest::test_ideal_3d()
00225 {
00226     MsqError err;
00227     MsqMatrix< 3, 3 > J;
00228     tet.MappingFunction3D::ideal( Sample( 3, 0 ), J, err );
00229     ASSERT_NO_ERROR( err );
00230     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( J ), 1e-6 );
00231 
00232     const Vector3D* verts = unit_edge_element( TETRAHEDRON );
00233     CPPUNIT_ASSERT( verts );
00234 
00235     size_t indices[4], num;
00236     MsqVector< 3 > coeff[4];
00237     tet.derivatives( Sample( 2, 0 ), NodeSet(), indices, coeff, num, err );
00238     ASSERT_NO_ERROR( err );
00239 
00240     MsqMatrix< 3, 3 > J_exp( 0.0 );  // zero
00241     for( size_t j = 0; j < num; ++j )
00242         J_exp += MsqVector< 3 >( verts[j].to_array() ) * transpose( coeff[j] );
00243     J_exp /= MBMesquite::cbrt( det( J_exp ) );
00244 
00245     // Matrices should be a rotation of each other.
00246     // First, calculate tentative rotation matrix
00247     MsqMatrix< 3, 3 > R = inverse( J_exp ) * J;
00248     // next check that it is a rotation
00249     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
00250     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
00251 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines