MOAB: Mesh Oriented datABase  (version 5.3.0)
TargetCalculatorTest.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) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file TargetCalculatorTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "UnitUtil.hpp"
00034 #include "TargetCalculator.hpp"
00035 #include "MsqError.hpp"
00036 #include "PatchData.hpp"
00037 #include "IdealElements.hpp"
00038 #include "LinearHexahedron.hpp"
00039 #include "LinearTriangle.hpp"
00040 #include "ReferenceMesh.hpp"
00041 #include <iostream>
00042 
00043 using namespace MBMesquite;
00044 
00045 double EPS    = 0.05;
00046 double EPSBIG = 0.5;
00047 
00048 class TargetCalculatorTest : public CppUnit::TestFixture
00049 {
00050   private:
00051     CPPUNIT_TEST_SUITE( TargetCalculatorTest );
00052 
00053     CPPUNIT_TEST( test_factor_2D );
00054     CPPUNIT_TEST( test_factor_surface );
00055     CPPUNIT_TEST( test_factor_3D );
00056     CPPUNIT_TEST( test_factor_2D_zero );
00057     CPPUNIT_TEST( test_factor_surface_zero );
00058     CPPUNIT_TEST( test_factor_3D_zero );
00059 
00060     CPPUNIT_TEST( test_size_2D );
00061     CPPUNIT_TEST( test_size_surface );
00062     CPPUNIT_TEST( test_size_3D );
00063 
00064     CPPUNIT_TEST( test_skew_2D );
00065     CPPUNIT_TEST( test_skew_surface );
00066     CPPUNIT_TEST( test_skew_3D );
00067 
00068     CPPUNIT_TEST( test_aspect_2D );
00069     CPPUNIT_TEST( test_aspect_surface );
00070     CPPUNIT_TEST( test_aspect_3D );
00071 
00072     CPPUNIT_TEST( test_shape_2D );
00073     CPPUNIT_TEST( test_shape_surface );
00074     CPPUNIT_TEST( test_shape_3D );
00075 
00076     CPPUNIT_TEST( test_ideal_skew_tri );
00077     CPPUNIT_TEST( test_ideal_skew_quad );
00078     CPPUNIT_TEST( test_ideal_skew_tet );
00079     CPPUNIT_TEST( test_ideal_skew_prism );
00080     CPPUNIT_TEST( test_ideal_skew_pyramid );
00081     CPPUNIT_TEST( test_ideal_skew_hex );
00082 
00083     CPPUNIT_TEST( test_ideal_shape_tri );
00084     CPPUNIT_TEST( test_ideal_shape_quad );
00085     CPPUNIT_TEST( test_ideal_shape_tet );
00086     CPPUNIT_TEST( test_ideal_shape_prism );
00087     CPPUNIT_TEST( test_ideal_shape_pyramid );
00088     CPPUNIT_TEST( test_ideal_shape_hex );
00089 
00090     CPPUNIT_TEST( test_new_orientatin_3D );
00091     CPPUNIT_TEST( test_new_orientatin_2D );
00092 
00093     CPPUNIT_TEST( test_new_aspect_3D );
00094     CPPUNIT_TEST( test_new_aspect_2D );
00095 
00096     CPPUNIT_TEST( test_jacobian_3D );
00097     CPPUNIT_TEST( test_jacobian_2D );
00098 
00099     CPPUNIT_TEST( test_get_refmesh_Jacobian_3D );
00100     CPPUNIT_TEST( test_get_refmesh_Jacobian_2D );
00101 
00102     CPPUNIT_TEST_SUITE_END();
00103 
00104     // define some matrices for use in testing.
00105     MsqMatrix< 3, 3 > V3D_Z45, V3D_X90, Q3D_45, D3D_123;
00106     MsqMatrix< 3, 2 > V2D_Z45, V2D_X90;
00107     MsqMatrix< 2, 2 > Q2D_45, D2D_21;
00108 
00109   public:
00110     TargetCalculatorTest();
00111     void setUp();
00112 
00113     void test_factor_2D();
00114     void test_factor_surface();
00115     void test_factor_3D();
00116     void test_factor_2D_zero();
00117     void test_factor_surface_zero();
00118     void test_factor_3D_zero();
00119 
00120     void test_size_2D();
00121     void test_size_surface();
00122     void test_size_3D();
00123 
00124     void test_skew_2D();
00125     void test_skew_surface();
00126     void test_skew_3D();
00127 
00128     void test_aspect_2D();
00129     void test_aspect_surface();
00130     void test_aspect_3D();
00131 
00132     void test_shape_2D();
00133     void test_shape_surface();
00134     void test_shape_3D();
00135 
00136     void test_ideal_skew_tri();
00137     void test_ideal_skew_quad();
00138     void test_ideal_skew_tet();
00139     void test_ideal_skew_prism();
00140     void test_ideal_skew_pyramid();
00141     void test_ideal_skew_hex();
00142 
00143     void test_ideal_shape_tri();
00144     void test_ideal_shape_quad();
00145     void test_ideal_shape_tet();
00146     void test_ideal_shape_prism();
00147     void test_ideal_shape_pyramid();
00148     void test_ideal_shape_hex();
00149 
00150     void test_new_orientatin_3D();
00151     void test_new_orientatin_2D();
00152 
00153     void test_new_aspect_3D();
00154     void test_new_aspect_2D();
00155 
00156     void test_jacobian_3D();
00157     void test_jacobian_2D();
00158 
00159     void test_get_refmesh_Jacobian_3D();
00160     void test_get_refmesh_Jacobian_2D();
00161 
00162     template < unsigned R, unsigned C >
00163     static inline void check_valid_V( MsqMatrix< R, C > V );
00164 
00165     template < unsigned D >
00166     static inline void check_valid_Q( MsqMatrix< D, D > Q );
00167 
00168     template < unsigned D >
00169     static inline void check_valid_delta( MsqMatrix< D, D > delta );
00170 };
00171 
00172 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TargetCalculatorTest, "TargetCalculatorTest" );
00173 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( TargetCalculatorTest, "Unit" );
00174 
00175 TargetCalculatorTest::TargetCalculatorTest()
00176 {
00177     const double cos45              = MSQ_SQRT_TWO / 2.0;
00178     const double rotation_3D_Z45[9] = { cos45, -cos45, 0, cos45, cos45, 0, 0, 0, 1 };
00179     const double rotation_2D_Z45[6] = { cos45, -cos45, cos45, cos45, 0, 0 };
00180 
00181     const double rotation_3D_X90[9] = { 1, 0, 0, 0, 0, -1, 0, 1, 0 };
00182     const double rotation_2D_X90[6] = { 1, 0, 0, 0, 0, 1 };
00183 
00184     const double rc45          = sqrt( cos45 );
00185     const double skew_2D_45[4] = { 1 / rc45, rc45, 0, rc45 };
00186     const double skew_3D_45[9] = { 1, cos45, cos45, 0, cos45, 1 - cos45, 0, 0, sqrt( MSQ_SQRT_TWO - 1 ) };
00187 
00188     const double aspect_2D_2x[4]  = { MSQ_SQRT_TWO, 0, 0, MSQ_SQRT_TWO / 2 };
00189     const double r6               = MBMesquite::cbrt( 1.0 / 6.0 );
00190     const double aspect_3D_123[9] = { r6, 0, 0, 0, 2 * r6, 0, 0, 0, 3 * r6 };
00191 
00192     V3D_Z45 = MsqMatrix< 3, 3 >( rotation_3D_Z45 );
00193     V3D_X90 = MsqMatrix< 3, 3 >( rotation_3D_X90 );
00194     Q3D_45  = MsqMatrix< 3, 3 >( skew_3D_45 );
00195     Q3D_45 *= 1 / Mesquite::cbrt( det( Q3D_45 ) );
00196     D3D_123 = MsqMatrix< 3, 3 >( aspect_3D_123 );
00197 
00198     V2D_Z45 = MsqMatrix< 3, 2 >( rotation_2D_Z45 );
00199     V2D_X90 = MsqMatrix< 3, 2 >( rotation_2D_X90 );
00200     Q2D_45  = MsqMatrix< 2, 2 >( skew_2D_45 );
00201     D2D_21  = MsqMatrix< 2, 2 >( aspect_2D_2x );
00202 }
00203 
00204 template < unsigned R, unsigned C >
00205 inline void TargetCalculatorTest::check_valid_V( MsqMatrix< R, C > V )
00206 {
00207     // check that it is a rotation
00208     MsqMatrix< C, C > I( 1.0 );
00209     ASSERT_MATRICES_EQUAL( I, transpose( V ) * V, EPS );
00210 }
00211 
00212 template < unsigned D >
00213 inline void TargetCalculatorTest::check_valid_Q( MsqMatrix< D, D > Q )
00214 {
00215     // must have unit determinant
00216     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( Q ), EPS );
00217 
00218     // must be upper-triangular
00219     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, Q( 1, 0 ), EPS );
00220     if( D == 3 )
00221     {
00222         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, Q( 2, 0 ), EPS );
00223         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, Q( 2, 1 ), EPS );
00224     }
00225 
00226     // columns must be of equal length
00227     CPPUNIT_ASSERT_DOUBLES_EQUAL( length( Q.column( 0 ) ), length( Q.column( 1 ) ), EPSBIG );
00228     if( D == 3 ) { CPPUNIT_ASSERT_DOUBLES_EQUAL( length( Q.column( 0 ) ), length( Q.column( 2 ) ), EPSBIG ); }
00229 
00230     // diagonal elements must be greater than zero
00231     CPPUNIT_ASSERT( Q( 0, 0 ) - EPS >= 0.0 );
00232     CPPUNIT_ASSERT( Q( 1, 1 ) - EPS >= 0.0 );
00233     if( D == 3 ) { CPPUNIT_ASSERT( Q( 2, 2 ) - EPS >= 0.0 ); }
00234 }
00235 
00236 template < unsigned D >
00237 inline void TargetCalculatorTest::check_valid_delta( MsqMatrix< D, D > delta )
00238 {
00239     // must have unit determinant
00240     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( delta ), EPS );
00241 
00242     // must be diagonal
00243     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, delta( 1, 0 ), EPS );
00244     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, delta( 0, 1 ), EPS );
00245     if( D == 3 )
00246     {
00247         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, delta( 2, 0 ), EPS );
00248         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, delta( 0, 2 ), EPS );
00249         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, delta( 1, 2 ), EPS );
00250         CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, delta( 2, 1 ), EPS );
00251     }
00252 }
00253 
00254 void TargetCalculatorTest::setUp()
00255 {
00256     // test the test: make sure we're testing with
00257     // valid matrix factors.
00258 
00259     check_valid_V( V3D_Z45 );
00260     check_valid_V( V3D_X90 );
00261     check_valid_V( V2D_Z45 );
00262     check_valid_V( V2D_X90 );
00263 
00264     check_valid_Q( Q3D_45 );
00265     check_valid_Q( Q2D_45 );
00266 
00267     check_valid_delta( D3D_123 );
00268     check_valid_delta( D2D_21 );
00269 }
00270 
00271 void TargetCalculatorTest::test_factor_2D()
00272 {
00273     MsqPrintError err( std::cout );
00274     MsqMatrix< 2, 2 > I( 1.0 ), W, V;
00275     double lambda;
00276     MsqMatrix< 2, 2 > Q, delta;
00277     bool valid;
00278 
00279     // first test with I
00280     valid = TargetCalculator::factor_2D( I, lambda, V, Q, delta, err );
00281     CPPUNIT_ASSERT( valid && !err );
00282     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, lambda, EPS );
00283     check_valid_V( V );
00284     check_valid_Q( Q );
00285     check_valid_delta( delta );
00286     W = lambda * V * Q * delta;
00287     ASSERT_MATRICES_EQUAL( I, W, EPS );
00288 
00289     // now test with 2*I
00290     valid = TargetCalculator::factor_2D( 2 * I, lambda, V, Q, delta, err );
00291     CPPUNIT_ASSERT( valid && !err );
00292     CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, lambda, EPS );
00293     check_valid_V( V );
00294     check_valid_Q( Q );
00295     check_valid_delta( delta );
00296     W = lambda * V * Q * delta;
00297     ASSERT_MATRICES_EQUAL( 2 * I, W, EPS );
00298 
00299     // now rotate the matrix about the X axis by 90 degrees
00300     MsqMatrix< 2, 2 > I_rot( 0.0 );
00301     I_rot( 0, 1 ) = 1.0;
00302     I_rot( 1, 0 ) = -1.0;
00303     valid         = TargetCalculator::factor_2D( I_rot, lambda, V, Q, delta, err );
00304     CPPUNIT_ASSERT( valid && !err );
00305     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, lambda, EPS );
00306     check_valid_V( V );
00307     check_valid_Q( Q );
00308     check_valid_delta( delta );
00309     W = lambda * V * Q * delta;
00310     ASSERT_MATRICES_EQUAL( I_rot, W, EPS );
00311 
00312     // change the aspect ratio
00313     MsqMatrix< 2, 2 > A( 1.0 );
00314     A( 1, 1 ) = 1.5;
00315     valid     = TargetCalculator::factor_2D( A, lambda, V, Q, delta, err );
00316     CPPUNIT_ASSERT( valid && !err );
00317     CPPUNIT_ASSERT_DOUBLES_EQUAL( pow( det( transpose( A ) * A ), 0.25 ), lambda, EPS );
00318     check_valid_V( V );
00319     check_valid_Q( Q );
00320     check_valid_delta( delta );
00321     W = lambda * V * Q * delta;
00322     ASSERT_MATRICES_EQUAL( A, W, EPS );
00323 
00324     // try an arbitrary matrix
00325     double w[4] = { 3, 1, 4, 2 };
00326     MsqMatrix< 2, 2 > W2( w );
00327     valid = TargetCalculator::factor_2D( W2, lambda, V, Q, delta, err );
00328     CPPUNIT_ASSERT( valid && !err );
00329     CPPUNIT_ASSERT_DOUBLES_EQUAL( pow( det( transpose( W2 ) * W2 ), 0.25 ), lambda, EPS );
00330     check_valid_V( V );
00331     check_valid_Q( Q );
00332     check_valid_delta( delta );
00333     W = lambda * V * Q * delta;
00334     ASSERT_MATRICES_EQUAL( W2, W, EPS );
00335 
00336     // try a more elaborate test
00337     const double e = exp( (double)1 );
00338     MsqMatrix< 2, 2 > Z45( V2D_Z45.data() );  // copy first two rows
00339     W2    = e * Z45 * Q2D_45 * D2D_21;
00340     valid = TargetCalculator::factor_2D( W2, lambda, V, Q, delta, err );
00341     CPPUNIT_ASSERT( valid && !err );
00342     CPPUNIT_ASSERT_DOUBLES_EQUAL( e, lambda, EPS );
00343     ASSERT_MATRICES_EQUAL( Z45, V, EPS );
00344     ASSERT_MATRICES_EQUAL( Q2D_45, Q, EPS );
00345     ASSERT_MATRICES_EQUAL( D2D_21, delta, EPS );
00346 }
00347 
00348 void TargetCalculatorTest::test_factor_surface()
00349 {
00350     MsqPrintError err( std::cout );
00351     MsqMatrix< 3, 2 > I( 1.0 ), W, V;
00352     double lambda;
00353     MsqMatrix< 2, 2 > Q, delta;
00354     bool valid;
00355 
00356     // first test with I
00357     valid = TargetCalculator::factor_surface( I, lambda, V, Q, delta, err );
00358     CPPUNIT_ASSERT( valid && !err );
00359     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, lambda, EPS );
00360     check_valid_V( V );
00361     check_valid_Q( Q );
00362     check_valid_delta( delta );
00363     W = lambda * V * Q * delta;
00364     ASSERT_MATRICES_EQUAL( I, W, EPS );
00365 
00366     // now test with 2*I
00367     valid = TargetCalculator::factor_surface( 2 * I, lambda, V, Q, delta, err );
00368     CPPUNIT_ASSERT( valid && !err );
00369     CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, lambda, EPS );
00370     check_valid_V( V );
00371     check_valid_Q( Q );
00372     check_valid_delta( delta );
00373     W = lambda * V * Q * delta;
00374     ASSERT_MATRICES_EQUAL( 2 * I, W, EPS );
00375 
00376     // now rotate the matrix about the X axis by 90 degrees
00377     MsqMatrix< 3, 2 > I_rot( 1.0 );
00378     I_rot( 1, 1 ) = 0.0;
00379     I_rot( 2, 1 ) = 1.0;
00380     valid         = TargetCalculator::factor_surface( I_rot, lambda, V, Q, delta, err );
00381     CPPUNIT_ASSERT( valid && !err );
00382     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, lambda, EPS );
00383     check_valid_V( V );
00384     check_valid_Q( Q );
00385     check_valid_delta( delta );
00386     W = lambda * V * Q * delta;
00387     ASSERT_MATRICES_EQUAL( I_rot, W, EPS );
00388 
00389     // change the aspect ratio
00390     MsqMatrix< 3, 2 > A( 1.0 );
00391     A( 1, 1 ) = 1.5;
00392     valid     = TargetCalculator::factor_surface( A, lambda, V, Q, delta, err );
00393     CPPUNIT_ASSERT( valid && !err );
00394     CPPUNIT_ASSERT_DOUBLES_EQUAL( pow( det( transpose( A ) * A ), 0.25 ), lambda, EPS );
00395     check_valid_V( V );
00396     check_valid_Q( Q );
00397     check_valid_delta( delta );
00398     W = lambda * V * Q * delta;
00399     ASSERT_MATRICES_EQUAL( A, W, EPS );
00400 
00401     // try an arbitrary matrix
00402     double w[6] = { 3, 1, 4, 2, -1, 5 };
00403     MsqMatrix< 3, 2 > W2( w );
00404     valid = TargetCalculator::factor_surface( W2, lambda, V, Q, delta, err );
00405     CPPUNIT_ASSERT( valid && !err );
00406     CPPUNIT_ASSERT_DOUBLES_EQUAL( pow( det( transpose( W2 ) * W2 ), 0.25 ), lambda, EPS );
00407     check_valid_V( V );
00408     check_valid_Q( Q );
00409     check_valid_delta( delta );
00410     W = lambda * V * Q * delta;
00411     ASSERT_MATRICES_EQUAL( W2, W, EPS );
00412 
00413     // try a more elaborate test
00414     const double e = exp( (double)1 );
00415     W2             = e * V2D_Z45 * Q2D_45 * D2D_21;
00416     valid          = TargetCalculator::factor_surface( W2, lambda, V, Q, delta, err );
00417     CPPUNIT_ASSERT( valid && !err );
00418     CPPUNIT_ASSERT_DOUBLES_EQUAL( e, lambda, EPS );
00419     ASSERT_MATRICES_EQUAL( V2D_Z45, V, EPS );
00420     ASSERT_MATRICES_EQUAL( Q2D_45, Q, EPS );
00421     ASSERT_MATRICES_EQUAL( D2D_21, delta, EPS );
00422 }
00423 
00424 void TargetCalculatorTest::test_factor_3D()
00425 {
00426     MsqPrintError err( std::cout );
00427     MsqMatrix< 3, 3 > I( 1.0 ), W, V, Q, delta;
00428     double lambda;
00429     bool valid;
00430 
00431     // first test with I
00432     valid = TargetCalculator::factor_3D( I, lambda, V, Q, delta, err );
00433     CPPUNIT_ASSERT( valid && !err );
00434     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, lambda, EPS );
00435     check_valid_V( V );
00436     check_valid_Q( Q );
00437     check_valid_delta( delta );
00438     W = lambda * V * Q * delta;
00439     ASSERT_MATRICES_EQUAL( I, W, EPS );
00440 
00441     // now test with 2*I
00442     valid = TargetCalculator::factor_3D( 2 * I, lambda, V, Q, delta, err );
00443     CPPUNIT_ASSERT( valid && !err );
00444     CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, lambda, EPS );
00445     check_valid_V( V );
00446     check_valid_Q( Q );
00447     check_valid_delta( delta );
00448     W = lambda * V * Q * delta;
00449     ASSERT_MATRICES_EQUAL( 2 * I, W, EPS );
00450 
00451     // now rotate the matrix about the X axis by 90 degrees
00452     MsqMatrix< 3, 3 > I_rot( 1.0 );
00453     I_rot( 1, 1 ) = 0.0;
00454     I_rot( 2, 1 ) = 1.0;
00455     I_rot( 1, 2 ) = -1.0;
00456     I_rot( 2, 2 ) = 0.0;
00457     valid         = TargetCalculator::factor_3D( I_rot, lambda, V, Q, delta, err );
00458     CPPUNIT_ASSERT( valid && !err );
00459     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, lambda, EPS );
00460     check_valid_V( V );
00461     check_valid_Q( Q );
00462     check_valid_delta( delta );
00463     W = lambda * V * Q * delta;
00464     ASSERT_MATRICES_EQUAL( I_rot, W, EPS );
00465 
00466     // change the aspect ratio
00467     MsqMatrix< 3, 3 > A( 1.0 );
00468     A( 1, 1 ) = 1.5;
00469     valid     = TargetCalculator::factor_3D( A, lambda, V, Q, delta, err );
00470     CPPUNIT_ASSERT( valid && !err );
00471     CPPUNIT_ASSERT_DOUBLES_EQUAL( MBMesquite::cbrt( det( A ) ), lambda, EPS );
00472     check_valid_V( V );
00473     check_valid_Q( Q );
00474     check_valid_delta( delta );
00475     W = lambda * V * Q * delta;
00476     ASSERT_MATRICES_EQUAL( A, W, EPS );
00477 
00478     // try an arbitrary matrix
00479     double vals[] = { 9, 8, 7, 1, 5, 4, 3, 2, 6 };
00480     MsqMatrix< 3, 3 > W2( vals );
00481     valid = TargetCalculator::factor_3D( W2, lambda, V, Q, delta, err );
00482     CPPUNIT_ASSERT( valid && !err );
00483     CPPUNIT_ASSERT_DOUBLES_EQUAL( MBMesquite::cbrt( det( W2 ) ), lambda, EPS );
00484     check_valid_V( V );
00485     check_valid_Q( Q );
00486     check_valid_delta( delta );
00487     W = lambda * V * Q * delta;
00488     ASSERT_MATRICES_EQUAL( W2, W, EPSBIG );
00489 
00490     // try a more elaborate test
00491     const double e = exp( (double)1 );
00492     W2             = e * V3D_Z45 * Q3D_45 * D3D_123;
00493     valid          = TargetCalculator::factor_3D( W2, lambda, V, Q, delta, err );
00494     CPPUNIT_ASSERT( valid && !err );
00495     CPPUNIT_ASSERT_DOUBLES_EQUAL( e, lambda, EPS );
00496     ASSERT_MATRICES_EQUAL( V3D_Z45, V, EPS );
00497     ASSERT_MATRICES_EQUAL( Q3D_45, Q, EPS );
00498     ASSERT_MATRICES_EQUAL( D3D_123, delta, EPS );
00499 }
00500 
00501 void TargetCalculatorTest::test_factor_2D_zero()
00502 {
00503     MsqPrintError err( std::cout );
00504     double lambda;
00505     MsqMatrix< 2, 2 > Q, delta;
00506     MsqMatrix< 2, 2 > V, Z( 0.0 );
00507     Z( 0, 0 ) = 1.0;
00508 
00509     bool valid = TargetCalculator::factor_2D( Z, lambda, V, Q, delta, err );
00510     ASSERT_NO_ERROR( err );
00511     CPPUNIT_ASSERT( !valid );
00512 }
00513 
00514 void TargetCalculatorTest::test_factor_surface_zero()
00515 {
00516     MsqPrintError err( std::cout );
00517     double lambda;
00518     MsqMatrix< 2, 2 > Q, delta;
00519     MsqMatrix< 3, 2 > V, Z( 0.0 );
00520     Z( 0, 0 ) = 1.0;
00521 
00522     bool valid = TargetCalculator::factor_surface( Z, lambda, V, Q, delta, err );
00523     ASSERT_NO_ERROR( err );
00524     CPPUNIT_ASSERT( !valid );
00525 }
00526 
00527 void TargetCalculatorTest::test_factor_3D_zero()
00528 {
00529     MsqPrintError err( std::cout );
00530     double lambda;
00531     MsqMatrix< 3, 3 > V, Q, delta, Z( 0.0 );
00532     Z( 0, 0 ) = 1.0;
00533 
00534     bool valid = TargetCalculator::factor_3D( Z, lambda, V, Q, delta, err );
00535     ASSERT_NO_ERROR( err );
00536     CPPUNIT_ASSERT( !valid );
00537 }
00538 
00539 void TargetCalculatorTest::test_size_2D()
00540 {
00541     MsqPrintError err( std::cout );
00542     double lambda;
00543     MsqMatrix< 2, 2 > V, Q, delta;
00544 
00545     MsqMatrix< 2, 2 > I( 1.0 );
00546     TargetCalculator::factor_2D( I, lambda, V, Q, delta, err );
00547     ASSERT_NO_ERROR( err );
00548     CPPUNIT_ASSERT_DOUBLES_EQUAL( lambda, TargetCalculator::size( I ), EPS );
00549 
00550     double w[] = { 3, 1, 4, 2 };
00551     MsqMatrix< 2, 2 > W( w );
00552     TargetCalculator::factor_2D( W, lambda, V, Q, delta, err );
00553     ASSERT_NO_ERROR( err );
00554     CPPUNIT_ASSERT_DOUBLES_EQUAL( lambda, TargetCalculator::size( W ), EPS );
00555 }
00556 
00557 void TargetCalculatorTest::test_size_surface()
00558 {
00559     MsqPrintError err( std::cout );
00560     double lambda;
00561     MsqMatrix< 3, 2 > V;
00562     MsqMatrix< 2, 2 > Q, delta;
00563 
00564     MsqMatrix< 3, 2 > I( 1.0 );
00565     TargetCalculator::factor_surface( I, lambda, V, Q, delta, err );
00566     ASSERT_NO_ERROR( err );
00567     CPPUNIT_ASSERT_DOUBLES_EQUAL( lambda, TargetCalculator::size( I ), EPS );
00568 
00569     double w[] = { 3, 1, 4, 2, -1, 5 };
00570     MsqMatrix< 3, 2 > W( w );
00571     TargetCalculator::factor_surface( W, lambda, V, Q, delta, err );
00572     ASSERT_NO_ERROR( err );
00573     CPPUNIT_ASSERT_DOUBLES_EQUAL( lambda, TargetCalculator::size( W ), EPS );
00574 }
00575 
00576 void TargetCalculatorTest::test_size_3D()
00577 {
00578     MsqPrintError err( std::cout );
00579     double lambda;
00580     MsqMatrix< 3, 3 > V, Q, delta;
00581 
00582     MsqMatrix< 3, 3 > I( 1.0 );
00583     TargetCalculator::factor_3D( I, lambda, V, Q, delta, err );
00584     ASSERT_NO_ERROR( err );
00585     CPPUNIT_ASSERT_DOUBLES_EQUAL( lambda, TargetCalculator::size( I ), EPS );
00586 
00587     double w[] = { 9, 8, 7, 1, 5, 4, 3, 2, 6 };
00588     MsqMatrix< 3, 3 > W( w );
00589     TargetCalculator::factor_3D( W, lambda, V, Q, delta, err );
00590     ASSERT_NO_ERROR( err );
00591     CPPUNIT_ASSERT_DOUBLES_EQUAL( lambda, TargetCalculator::size( W ), EPS );
00592 }
00593 
00594 void TargetCalculatorTest::test_skew_2D()
00595 {
00596     MsqPrintError err( std::cout );
00597     double lambda;
00598     MsqMatrix< 2, 2 > V, Q, delta;
00599 
00600     MsqMatrix< 2, 2 > I( 1.0 );
00601     TargetCalculator::factor_2D( I, lambda, V, Q, delta, err );
00602     ASSERT_NO_ERROR( err );
00603     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( I ), EPS );
00604 
00605     double r[] = { 0, -1, 2, 0 };
00606     MsqMatrix< 2, 2 > R( r );
00607     TargetCalculator::factor_2D( R, lambda, V, Q, delta, err );
00608     ASSERT_NO_ERROR( err );
00609     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( R ), EPS );
00610 
00611     double w[] = { 3, 1, 4, 2 };
00612     MsqMatrix< 2, 2 > W( w );
00613     TargetCalculator::factor_2D( W, lambda, V, Q, delta, err );
00614     ASSERT_NO_ERROR( err );
00615     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( W ), EPS );
00616 
00617     W = MsqMatrix< 2, 2 >( 6 ) * Q2D_45;
00618     ASSERT_MATRICES_EQUAL( Q2D_45, TargetCalculator::skew( W ), EPS );
00619 }
00620 
00621 void TargetCalculatorTest::test_skew_surface()
00622 {
00623     MsqPrintError err( std::cout );
00624     double lambda;
00625     MsqMatrix< 3, 2 > V;
00626     MsqMatrix< 2, 2 > Q, delta;
00627 
00628     MsqMatrix< 3, 2 > I( 1.0 );
00629     TargetCalculator::factor_surface( I, lambda, V, Q, delta, err );
00630     ASSERT_NO_ERROR( err );
00631     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( I ), EPS );
00632 
00633     double r[] = { 1, 0, 0, 0, 0, 2 };
00634     MsqMatrix< 3, 2 > R( r );
00635     TargetCalculator::factor_surface( R, lambda, V, Q, delta, err );
00636     ASSERT_NO_ERROR( err );
00637     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( R ), EPS );
00638 
00639     double w[] = { 3, 1, 4, 2, -1, 5 };
00640     MsqMatrix< 3, 2 > W( w );
00641     TargetCalculator::factor_surface( W, lambda, V, Q, delta, err );
00642     ASSERT_NO_ERROR( err );
00643     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( W ), EPS );
00644 
00645     W = MsqMatrix< 3, 2 >( 6 ) * Q2D_45;
00646     ASSERT_MATRICES_EQUAL( Q2D_45, TargetCalculator::skew( W ), EPS );
00647 }
00648 
00649 void TargetCalculatorTest::test_skew_3D()
00650 {
00651     MsqPrintError err( std::cout );
00652     double lambda;
00653     MsqMatrix< 3, 3 > V, Q, delta;
00654 
00655     MsqMatrix< 3, 3 > I( 1.0 );
00656     TargetCalculator::factor_3D( I, lambda, V, Q, delta, err );
00657     ASSERT_NO_ERROR( err );
00658     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( I ), EPS );
00659 
00660     double r[] = { 0.8, 0.0, 0.0, 0.0, 0.0, 1.13, 0.0, -0.5, 0.0 };
00661     MsqMatrix< 3, 3 > R( r );
00662     TargetCalculator::factor_3D( R, lambda, V, Q, delta, err );
00663     ASSERT_NO_ERROR( err );
00664     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( R ), EPS );
00665 
00666     double w[] = { 9, 8, 7, 1, 5, 4, 3, 2, 6 };
00667     MsqMatrix< 3, 3 > W( w );
00668     TargetCalculator::factor_3D( W, lambda, V, Q, delta, err );
00669     ASSERT_NO_ERROR( err );
00670     ASSERT_MATRICES_EQUAL( Q, TargetCalculator::skew( W ), 2 * EPS );
00671 
00672     W = MsqMatrix< 3, 3 >( 2.1 ) * Q3D_45;
00673     ASSERT_MATRICES_EQUAL( Q3D_45, TargetCalculator::skew( W ), EPS );
00674 }
00675 
00676 void TargetCalculatorTest::test_aspect_2D()
00677 {
00678     MsqPrintError err( std::cout );
00679     double lambda;
00680     MsqMatrix< 2, 2 > V, Q, delta;
00681 
00682     MsqMatrix< 2, 2 > I( 1.0 );
00683     TargetCalculator::factor_2D( I, lambda, V, Q, delta, err );
00684     ASSERT_NO_ERROR( err );
00685     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( I ), EPS );
00686 
00687     double r[] = { 0, -1, 2, 0 };
00688     MsqMatrix< 2, 2 > R( r );
00689     TargetCalculator::factor_2D( R, lambda, V, Q, delta, err );
00690     ASSERT_NO_ERROR( err );
00691     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( R ), EPS );
00692 
00693     double w[] = { 3, 1, 4, 2 };
00694     MsqMatrix< 2, 2 > W( w );
00695     TargetCalculator::factor_2D( W, lambda, V, Q, delta, err );
00696     ASSERT_NO_ERROR( err );
00697     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( W ), EPS );
00698 
00699     W = MsqMatrix< 2, 2 >( 6 ) * Q2D_45;
00700     ASSERT_MATRICES_EQUAL( I, TargetCalculator::aspect( W ), EPS );
00701 }
00702 
00703 void TargetCalculatorTest::test_aspect_surface()
00704 {
00705     MsqPrintError err( std::cout );
00706     double lambda;
00707     MsqMatrix< 3, 2 > V;
00708     MsqMatrix< 2, 2 > Q, delta;
00709 
00710     MsqMatrix< 3, 2 > I( 1.0 );
00711     TargetCalculator::factor_surface( I, lambda, V, Q, delta, err );
00712     ASSERT_NO_ERROR( err );
00713     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( I ), EPS );
00714 
00715     double r[] = { 1, 0, 0, 0, 0, 2 };
00716     MsqMatrix< 3, 2 > R( r );
00717     TargetCalculator::factor_surface( R, lambda, V, Q, delta, err );
00718     ASSERT_NO_ERROR( err );
00719     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( R ), EPS );
00720 
00721     double w[] = { 3, 1, 4, 2, -1, 5 };
00722     MsqMatrix< 3, 2 > W( w );
00723     TargetCalculator::factor_surface( W, lambda, V, Q, delta, err );
00724     ASSERT_NO_ERROR( err );
00725     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( W ), EPS );
00726 
00727     W = MsqMatrix< 3, 2 >( 6 ) * Q2D_45;
00728     ASSERT_MATRICES_EQUAL( ( MsqMatrix< 2, 2 >( 1.0 ) ), TargetCalculator::aspect( W ), EPS );
00729 }
00730 
00731 void TargetCalculatorTest::test_aspect_3D()
00732 {
00733     MsqPrintError err( std::cout );
00734     double lambda;
00735     MsqMatrix< 3, 3 > V, Q, delta;
00736 
00737     MsqMatrix< 3, 3 > I( 1.0 );
00738     TargetCalculator::factor_3D( I, lambda, V, Q, delta, err );
00739     ASSERT_NO_ERROR( err );
00740     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( I ), EPS );
00741 
00742     double r[] = { 0.8, 0.0, 0.0, 0.0, 0.0, 1.13, 0.0, -0.5, 0.0 };
00743     MsqMatrix< 3, 3 > R( r );
00744     TargetCalculator::factor_3D( R, lambda, V, Q, delta, err );
00745     ASSERT_NO_ERROR( err );
00746     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( R ), EPS );
00747 
00748     double w[] = { 9, 8, 7, 1, 5, 4, 3, 2, 6 };
00749     MsqMatrix< 3, 3 > W( w );
00750     TargetCalculator::factor_3D( W, lambda, V, Q, delta, err );
00751     ASSERT_NO_ERROR( err );
00752     ASSERT_MATRICES_EQUAL( delta, TargetCalculator::aspect( W ), 2 * EPS );
00753 
00754     W = MsqMatrix< 3, 3 >( 2.1 ) * Q3D_45;
00755     ASSERT_MATRICES_EQUAL( I, TargetCalculator::aspect( W ), EPS );
00756 }
00757 
00758 void TargetCalculatorTest::test_shape_2D()
00759 {
00760     MsqPrintError err( std::cout );
00761     double lambda;
00762     MsqMatrix< 2, 2 > V, Q, delta;
00763 
00764     MsqMatrix< 2, 2 > I( 1.0 );
00765     TargetCalculator::factor_2D( I, lambda, V, Q, delta, err );
00766     ASSERT_NO_ERROR( err );
00767     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( I ), EPS );
00768 
00769     double r[] = { 0, -1, 2, 0 };
00770     MsqMatrix< 2, 2 > R( r );
00771     TargetCalculator::factor_2D( R, lambda, V, Q, delta, err );
00772     ASSERT_NO_ERROR( err );
00773     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( R ), EPS );
00774 
00775     double w[] = { 3, 1, 4, 2 };
00776     MsqMatrix< 2, 2 > W( w );
00777     TargetCalculator::factor_2D( W, lambda, V, Q, delta, err );
00778     ASSERT_NO_ERROR( err );
00779     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( W ), EPS );
00780 
00781     W = MsqMatrix< 2, 2 >( 6 ) * Q2D_45;
00782     ASSERT_MATRICES_EQUAL( Q2D_45, TargetCalculator::shape( W ), EPS );
00783 }
00784 
00785 void TargetCalculatorTest::test_shape_surface()
00786 {
00787     MsqPrintError err( std::cout );
00788     double lambda;
00789     MsqMatrix< 3, 2 > V;
00790     MsqMatrix< 2, 2 > Q, delta;
00791 
00792     MsqMatrix< 3, 2 > I( 1.0 );
00793     TargetCalculator::factor_surface( I, lambda, V, Q, delta, err );
00794     ASSERT_NO_ERROR( err );
00795     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( I ), EPS );
00796 
00797     double r[] = { 1, 0, 0, 0, 0, 2 };
00798     MsqMatrix< 3, 2 > R( r );
00799     TargetCalculator::factor_surface( R, lambda, V, Q, delta, err );
00800     ASSERT_NO_ERROR( err );
00801     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( R ), EPS );
00802 
00803     double w[] = { 3, 1, 4, 2, -1, 5 };
00804     MsqMatrix< 3, 2 > W( w );
00805     TargetCalculator::factor_surface( W, lambda, V, Q, delta, err );
00806     ASSERT_NO_ERROR( err );
00807     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( W ), EPS );
00808 
00809     W = MsqMatrix< 3, 2 >( 6 ) * Q2D_45;
00810     ASSERT_MATRICES_EQUAL( Q2D_45, TargetCalculator::shape( W ), EPS );
00811 }
00812 
00813 void TargetCalculatorTest::test_shape_3D()
00814 {
00815     MsqPrintError err( std::cout );
00816     double lambda;
00817     MsqMatrix< 3, 3 > V, Q, delta;
00818 
00819     MsqMatrix< 3, 3 > I( 1.0 );
00820     TargetCalculator::factor_3D( I, lambda, V, Q, delta, err );
00821     ASSERT_NO_ERROR( err );
00822     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( I ), EPS );
00823 
00824     double r[] = { 0.8, 0.0, 0.0, 0.0, 0.0, 1.13, 0.0, -0.5, 0.0 };
00825     MsqMatrix< 3, 3 > R( r );
00826     TargetCalculator::factor_3D( R, lambda, V, Q, delta, err );
00827     ASSERT_NO_ERROR( err );
00828     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( R ), EPS );
00829 
00830     double w[] = { 9, 8, 7, 1, 5, 4, 3, 2, 6 };
00831     MsqMatrix< 3, 3 > W( w );
00832     TargetCalculator::factor_3D( W, lambda, V, Q, delta, err );
00833     ASSERT_NO_ERROR( err );
00834     ASSERT_MATRICES_EQUAL( Q * delta, TargetCalculator::shape( W ), 2 * EPS );
00835 
00836     W = MsqMatrix< 3, 3 >( 2.1 ) * Q3D_45;
00837     ASSERT_MATRICES_EQUAL( Q3D_45, TargetCalculator::shape( W ), EPS );
00838 }
00839 
00840 void TargetCalculatorTest::test_ideal_skew_tri()
00841 {
00842     MsqError err;
00843     PatchData pd;
00844     MsqMatrix< 2, 2 > Q, Q2;
00845 
00846     TargetCalculator::ideal_skew_2D( TRIANGLE, Sample( 0, 0 ), pd, Q, err );
00847     ASSERT_NO_ERROR( err );
00848     check_valid_Q( Q );
00849 
00850     TargetCalculator::ideal_skew_2D( TRIANGLE, Sample( 2, 0 ), pd, Q2, err );
00851     ASSERT_NO_ERROR( err );
00852     ASSERT_MATRICES_EQUAL( Q, Q2, EPS );
00853 
00854     const Vector3D* coords = unit_edge_element( TRIANGLE );
00855     MsqMatrix< 3, 2 > W;
00856     TargetCalculator::jacobian_2D( pd, TRIANGLE, 3, Sample( 0, 0 ), coords, W, err );
00857     ASSERT_NO_ERROR( err );
00858     Q = TargetCalculator::skew( W );
00859     ASSERT_MATRICES_EQUAL( Q, Q2, EPS );
00860 }
00861 
00862 void TargetCalculatorTest::test_ideal_skew_quad()
00863 {
00864     MsqError err;
00865     PatchData pd;
00866     MsqMatrix< 2, 2 > Q;
00867 
00868     TargetCalculator::ideal_skew_2D( QUADRILATERAL, Sample( 0, 0 ), pd, Q, err );
00869     ASSERT_NO_ERROR( err );
00870     ASSERT_IDENTITY_MATRIX( Q );
00871 
00872     TargetCalculator::ideal_skew_2D( QUADRILATERAL, Sample( 2, 0 ), pd, Q, err );
00873     ASSERT_NO_ERROR( err );
00874     ASSERT_IDENTITY_MATRIX( Q );
00875 }
00876 
00877 void TargetCalculatorTest::test_ideal_skew_tet()
00878 {
00879     MsqError err;
00880     PatchData pd;
00881     MsqMatrix< 3, 3 > Q, Q2;
00882 
00883     TargetCalculator::ideal_skew_3D( TETRAHEDRON, Sample( 0, 0 ), pd, Q, err );
00884     ASSERT_NO_ERROR( err );
00885     check_valid_Q( Q );
00886 
00887     TargetCalculator::ideal_skew_3D( TETRAHEDRON, Sample( 3, 0 ), pd, Q2, err );
00888     ASSERT_NO_ERROR( err );
00889     ASSERT_MATRICES_EQUAL( Q, Q2, EPS );
00890 
00891     const Vector3D* coords = unit_edge_element( TETRAHEDRON );
00892     MsqMatrix< 3, 3 > W;
00893     TargetCalculator::jacobian_3D( pd, TETRAHEDRON, 4, Sample( 0, 0 ), coords, W, err );
00894     ASSERT_NO_ERROR( err );
00895     Q = TargetCalculator::skew( W );
00896     ASSERT_MATRICES_EQUAL( Q, Q2, EPS );
00897 }
00898 
00899 void TargetCalculatorTest::test_ideal_skew_prism()
00900 {
00901     Sample points[] = {
00902         Sample( 0, 0 ), Sample( 0, 1 ), Sample( 0, 2 ), Sample( 0, 3 ), Sample( 2, 0 ), Sample( 3, 0 )
00903     };
00904     const int num_pts = sizeof( points ) / sizeof( points[0] );
00905 
00906     MsqError err;
00907     PatchData pd;
00908     MsqMatrix< 3, 3 > Q, W;
00909     const Vector3D* coords = unit_edge_element( PRISM );
00910 
00911     for( int i = 0; i < num_pts; ++i )
00912     {
00913 
00914         TargetCalculator::ideal_skew_3D( PRISM, points[i], pd, Q, err );
00915         ASSERT_NO_ERROR( err );
00916         check_valid_Q( Q );
00917 
00918         TargetCalculator::jacobian_3D( pd, PRISM, 6, points[i], coords, W, err );
00919         ASSERT_NO_ERROR( err );
00920 
00921         ASSERT_MATRICES_EQUAL( TargetCalculator::skew( W ), Q, EPS );
00922     }
00923 }
00924 
00925 void TargetCalculatorTest::test_ideal_skew_pyramid()
00926 {
00927     Sample points[] = {
00928         Sample( 0, 0 ), Sample( 0, 1 ), Sample( 0, 2 ), Sample( 0, 3 ), Sample( 2, 0 ), Sample( 3, 0 )
00929     };
00930     const int num_pts = sizeof( points ) / sizeof( points[0] );
00931 
00932     MsqError err;
00933     PatchData pd;
00934     MsqMatrix< 3, 3 > Q, W;
00935     const Vector3D* coords = unit_edge_element( PYRAMID, true );
00936 
00937     for( int i = 0; i < num_pts; ++i )
00938     {
00939 
00940         TargetCalculator::ideal_skew_3D( PYRAMID, points[i], pd, Q, err );
00941         ASSERT_NO_ERROR( err );
00942         check_valid_Q( Q );
00943 
00944         TargetCalculator::jacobian_3D( pd, PYRAMID, 5, points[i], coords, W, err );
00945         ASSERT_NO_ERROR( err );
00946 
00947         ASSERT_MATRICES_EQUAL( TargetCalculator::skew( W ), Q, EPS );
00948     }
00949 }
00950 
00951 void TargetCalculatorTest::test_ideal_skew_hex()
00952 {
00953     MsqError err;
00954     PatchData pd;
00955     MsqMatrix< 3, 3 > Q;
00956 
00957     TargetCalculator::ideal_skew_3D( HEXAHEDRON, Sample( 0, 0 ), pd, Q, err );
00958     ASSERT_NO_ERROR( err );
00959     ASSERT_IDENTITY_MATRIX( Q );
00960 
00961     TargetCalculator::ideal_skew_3D( HEXAHEDRON, Sample( 3, 0 ), pd, Q, err );
00962     ASSERT_NO_ERROR( err );
00963     ASSERT_IDENTITY_MATRIX( Q );
00964 }
00965 
00966 void TargetCalculatorTest::test_ideal_shape_tri()
00967 {
00968     // Ideal triangle should have Aspect (i.e. delta) == identity,
00969     // so shape should be equal to skew.
00970 
00971     MsqError err;
00972     PatchData pd;
00973     MsqMatrix< 2, 2 > Q, Q2;
00974 
00975     TargetCalculator::ideal_shape_2D( TRIANGLE, Sample( 0, 0 ), pd, Q, err );
00976     ASSERT_NO_ERROR( err );
00977     TargetCalculator::ideal_skew_2D( TRIANGLE, Sample( 0, 0 ), pd, Q2, err );
00978     ASSERT_NO_ERROR( err );
00979     ASSERT_MATRICES_EQUAL( Q2, Q, EPS );
00980 
00981     TargetCalculator::ideal_shape_2D( TRIANGLE, Sample( 2, 0 ), pd, Q, err );
00982     ASSERT_NO_ERROR( err );
00983     TargetCalculator::ideal_skew_2D( TRIANGLE, Sample( 2, 0 ), pd, Q2, err );
00984     ASSERT_NO_ERROR( err );
00985     ASSERT_MATRICES_EQUAL( Q2, Q, EPS );
00986 }
00987 
00988 void TargetCalculatorTest::test_ideal_shape_quad()
00989 {
00990     MsqError err;
00991     PatchData pd;
00992     MsqMatrix< 2, 2 > Q;
00993 
00994     TargetCalculator::ideal_shape_2D( QUADRILATERAL, Sample( 0, 0 ), pd, Q, err );
00995     ASSERT_NO_ERROR( err );
00996     ASSERT_IDENTITY_MATRIX( Q );
00997 
00998     TargetCalculator::ideal_shape_2D( QUADRILATERAL, Sample( 2, 0 ), pd, Q, err );
00999     ASSERT_NO_ERROR( err );
01000     ASSERT_IDENTITY_MATRIX( Q );
01001 }
01002 
01003 void TargetCalculatorTest::test_ideal_shape_tet()
01004 {
01005     // Ideal tetrahedron should have Aspect (i.e. delta) == identity,
01006     // so shape should be equal to skew.
01007 
01008     MsqError err;
01009     PatchData pd;
01010     MsqMatrix< 3, 3 > Q, Q2;
01011 
01012     TargetCalculator::ideal_shape_3D( TETRAHEDRON, Sample( 0, 0 ), pd, Q, err );
01013     ASSERT_NO_ERROR( err );
01014     TargetCalculator::ideal_skew_3D( TETRAHEDRON, Sample( 0, 0 ), pd, Q2, err );
01015     ASSERT_NO_ERROR( err );
01016     ASSERT_MATRICES_EQUAL( Q2, Q, EPS );
01017 
01018     TargetCalculator::ideal_shape_3D( TETRAHEDRON, Sample( 3, 0 ), pd, Q, err );
01019     ASSERT_NO_ERROR( err );
01020     TargetCalculator::ideal_skew_3D( TETRAHEDRON, Sample( 3, 0 ), pd, Q2, err );
01021     ASSERT_NO_ERROR( err );
01022     ASSERT_MATRICES_EQUAL( Q2, Q, EPS );
01023 }
01024 
01025 void TargetCalculatorTest::test_ideal_shape_prism()
01026 {
01027     // Ideal wedge should have Aspect (i.e. delta) == identity,
01028     // so shape should be equal to skew.
01029 
01030     Sample points[] = {
01031         Sample( 0, 0 ), Sample( 0, 1 ), Sample( 0, 2 ), Sample( 0, 3 ), Sample( 2, 0 ), Sample( 3, 0 )
01032     };
01033     const int num_pts = sizeof( points ) / sizeof( points[0] );
01034 
01035     MsqError err;
01036     PatchData pd;
01037     MsqMatrix< 3, 3 > Q, W;
01038 
01039     for( int i = 0; i < num_pts; ++i )
01040     {
01041         TargetCalculator::ideal_shape_3D( PRISM, points[i], pd, Q, err );
01042         ASSERT_NO_ERROR( err );
01043         TargetCalculator::ideal_skew_3D( PRISM, points[i], pd, W, err );
01044         ASSERT_NO_ERROR( err );
01045         ASSERT_MATRICES_EQUAL( W, Q, EPS );
01046     }
01047 }
01048 
01049 void TargetCalculatorTest::test_ideal_shape_pyramid()
01050 {
01051     Sample points[] = {
01052         Sample( 0, 0 ), Sample( 0, 1 ), Sample( 0, 2 ), Sample( 0, 3 ), Sample( 2, 0 ), Sample( 3, 0 )
01053     };
01054     const int num_pts = sizeof( points ) / sizeof( points[0] );
01055 
01056     MsqError err;
01057     PatchData pd;
01058     MsqMatrix< 3, 3 > Q, W;
01059     const Vector3D* coords = unit_edge_element( PYRAMID, true );
01060 
01061     for( int i = 0; i < num_pts; ++i )
01062     {
01063 
01064         TargetCalculator::ideal_shape_3D( PYRAMID, points[i], pd, Q, err );
01065         ASSERT_NO_ERROR( err );
01066 
01067         TargetCalculator::jacobian_3D( pd, PYRAMID, 5, points[i], coords, W, err );
01068         ASSERT_NO_ERROR( err );
01069 
01070         ASSERT_MATRICES_EQUAL( TargetCalculator::shape( W ), Q, EPS );
01071     }
01072 }
01073 
01074 void TargetCalculatorTest::test_ideal_shape_hex()
01075 {
01076     MsqError err;
01077     PatchData pd;
01078     MsqMatrix< 3, 3 > Q;
01079 
01080     TargetCalculator::ideal_shape_3D( HEXAHEDRON, Sample( 0, 0 ), pd, Q, err );
01081     ASSERT_NO_ERROR( err );
01082     ASSERT_IDENTITY_MATRIX( Q );
01083 
01084     TargetCalculator::ideal_shape_3D( HEXAHEDRON, Sample( 3, 0 ), pd, Q, err );
01085     ASSERT_NO_ERROR( err );
01086     ASSERT_IDENTITY_MATRIX( Q );
01087 }
01088 
01089 void TargetCalculatorTest::test_new_orientatin_3D()
01090 {
01091     MsqVector< 3 > x( 0.0 ), y( 0.0 ), z( 0.0 );
01092     x[0] = y[1] = z[2] = 1;
01093 
01094     MsqMatrix< 3, 3 > V = TargetCalculator::new_orientation_3D( x, y );
01095     check_valid_V( V );
01096     ASSERT_MATRICES_EQUAL( z, V.column( 0 ) * V.column( 1 ), EPS );
01097     ASSERT_MATRICES_EQUAL( z, V.column( 2 ), EPS );
01098 }
01099 
01100 void TargetCalculatorTest::test_new_orientatin_2D()
01101 {
01102     MsqVector< 3 > x( 0.0 ), y( 0.0 ), z( 0.0 );
01103     x[0] = y[1] = z[2] = 1;
01104 
01105     MsqMatrix< 3, 2 > V = TargetCalculator::new_orientation_2D( x, y );
01106     check_valid_V( V );
01107     ASSERT_MATRICES_EQUAL( z, V.column( 0 ) * V.column( 1 ), EPS );
01108 }
01109 
01110 void TargetCalculatorTest::test_new_aspect_3D()
01111 {
01112     MsqVector< 3 > r;
01113     r[0]                    = 2;
01114     r[1]                    = 4;
01115     r[2]                    = 6;
01116     MsqMatrix< 3, 3 > delta = TargetCalculator::new_aspect_3D( r );
01117     check_valid_delta( delta );
01118     CPPUNIT_ASSERT_DOUBLES_EQUAL( r[0] / r[1], delta( 0, 0 ) / delta( 1, 1 ), EPS );
01119     CPPUNIT_ASSERT_DOUBLES_EQUAL( r[0] / r[2], delta( 0, 0 ) / delta( 2, 2 ), EPS );
01120 }
01121 
01122 void TargetCalculatorTest::test_new_aspect_2D()
01123 {
01124     MsqVector< 2 > r;
01125     r[0]                    = 7;
01126     r[1]                    = 3;
01127     MsqMatrix< 2, 2 > delta = TargetCalculator::new_aspect_2D( r );
01128     check_valid_delta( delta );
01129     CPPUNIT_ASSERT_DOUBLES_EQUAL( r[0] / r[1], delta( 0, 0 ) / delta( 1, 1 ), EPS );
01130 }
01131 
01132 void TargetCalculatorTest::test_jacobian_3D()
01133 {
01134     MsqError err;
01135     PatchData pd;
01136     LinearHexahedron map;
01137     const Vector3D hex_coords[] = {
01138         Vector3D( 0, 0, 0 ), Vector3D( 2, 0, 0 ), Vector3D( 3, 2, 0 ), Vector3D( 1, 2, 0 ),
01139         Vector3D( 1, 0, 1 ), Vector3D( 3, 0, 2 ), Vector3D( 3, 2, 2 ), Vector3D( 2, 2, 2 )
01140     };
01141     const MsqVector< 3 >* coords = reinterpret_cast< const MsqVector< 3 >* >( hex_coords );
01142 
01143     Sample pts[]      = { Sample( 0, 0 ), Sample( 1, 6 ), Sample( 3, 0 ) };
01144     const int num_pts = sizeof( pts ) / sizeof( pts[0] );
01145     for( int i = 0; i < num_pts; ++i )
01146     {
01147         MsqMatrix< 3, 3 > J( 0.0 ), W;
01148         size_t indices[8], n;
01149         MsqVector< 3 > derivs[8];
01150         map.derivatives( pts[i], NodeSet(), indices, derivs, n, err );
01151         ASSERT_NO_ERROR( err );
01152         for( size_t j = 0; j < n; ++j )
01153             J += outer( coords[indices[j]], derivs[j] );
01154 
01155         TargetCalculator::jacobian_3D( pd, HEXAHEDRON, 8, pts[i], hex_coords, W, err );
01156         ASSERT_NO_ERROR( err );
01157 
01158         ASSERT_MATRICES_EQUAL( J, W, EPS );
01159     }
01160 }
01161 
01162 void TargetCalculatorTest::test_jacobian_2D()
01163 {
01164     MsqError err;
01165     PatchData pd;
01166     LinearTriangle map;
01167     const Vector3D tri_coords[]  = { Vector3D( 0, 0, 0 ), Vector3D( 2, 0, 0 ), Vector3D( 0, 0, 1 ) };
01168     const MsqVector< 3 >* coords = reinterpret_cast< const MsqVector< 3 >* >( tri_coords );
01169 
01170     Sample pts[]      = { Sample( 0, 0 ), Sample( 1, 1 ), Sample( 2, 0 ) };
01171     const int num_pts = sizeof( pts ) / sizeof( pts[0] );
01172     for( int i = 0; i < num_pts; ++i )
01173     {
01174         MsqMatrix< 3, 2 > J( 0.0 ), W;
01175         size_t indices[3], n;
01176         MsqVector< 2 > derivs[3];
01177         map.derivatives( pts[i], NodeSet(), indices, derivs, n, err );
01178         ASSERT_NO_ERROR( err );
01179         for( size_t j = 0; j < n; ++j )
01180             J += outer( coords[indices[j]], derivs[j] );
01181 
01182         TargetCalculator::jacobian_2D( pd, TRIANGLE, 3, pts[i], tri_coords, W, err );
01183         ASSERT_NO_ERROR( err );
01184 
01185         ASSERT_MATRICES_EQUAL( J, W, EPS );
01186     }
01187 }
01188 
01189 class DummyRefMesh : public ReferenceMeshInterface
01190 {
01191     std::vector< Vector3D > mCoords;
01192 
01193   public:
01194     DummyRefMesh( const Vector3D* coords, int num_verts ) : mCoords( coords, coords + num_verts ) {}
01195 
01196     void get_reference_vertex_coordinates( const Mesh::VertexHandle* handles, const size_t num_vertices,
01197                                            Vector3D* coords, MsqError& err )
01198     {
01199         const size_t* indices = reinterpret_cast< const size_t* >( handles );
01200         for( size_t i = 0; i < num_vertices; ++i )
01201         {
01202             if( i >= mCoords.size() )
01203             {
01204                 MSQ_SETERR( err )( MsqError::INVALID_ARG );
01205                 return;
01206             }
01207             coords[i] = mCoords[indices[i]];
01208         }
01209     }
01210 };
01211 
01212 void TargetCalculatorTest::test_get_refmesh_Jacobian_3D()
01213 {
01214     MsqError err;
01215 
01216     const Vector3D hex_coords[] = {
01217         Vector3D( 0, 0, 0 ), Vector3D( 2, 0, 0 ), Vector3D( 3, 2, 0 ), Vector3D( 1, 2, 0 ),
01218         Vector3D( 1, 0, 1 ), Vector3D( 3, 0, 2 ), Vector3D( 3, 2, 2 ), Vector3D( 2, 2, 2 )
01219     };
01220     DummyRefMesh ref_mesh( hex_coords, 8 );
01221 
01222     const Vector3D rect_coords[] = {
01223         Vector3D( 0, 0, 0 ), Vector3D( 1, 0, 0 ), Vector3D( 1, 1, 0 ), Vector3D( 0, 1, 0 ),
01224         Vector3D( 0, 0, 5 ), Vector3D( 1, 0, 5 ), Vector3D( 1, 1, 5 ), Vector3D( 0, 1, 5 )
01225     };
01226     size_t conn[] = { 0, 1, 2, 3, 4, 5, 6, 7 };
01227     bool fixed[8];
01228     std::fill( fixed, fixed + sizeof( fixed ) / sizeof( fixed[0] ), false );
01229     PatchData pd;
01230     pd.fill( 8, rect_coords[0].to_array(), 1, HEXAHEDRON, conn, fixed, err );
01231     ASSERT_NO_ERROR( err );
01232 
01233     Sample pts[]      = { Sample( 0, 0 ), Sample( 1, 6 ), Sample( 3, 0 ) };
01234     const int num_pts = sizeof( pts ) / sizeof( pts[0] );
01235     for( int i = 0; i < num_pts; ++i )
01236     {
01237         MsqMatrix< 3, 3 > J, W;
01238         TargetCalculator::jacobian_3D( pd, HEXAHEDRON, 8, pts[i], hex_coords, J, err );
01239         ASSERT_NO_ERROR( err );
01240 
01241         TargetCalculator::get_refmesh_Jacobian_3D( &ref_mesh, pd, 0, pts[i], W, err );
01242         ASSERT_NO_ERROR( err );
01243 
01244         ASSERT_MATRICES_EQUAL( J, W, EPS );
01245     }
01246 }
01247 
01248 void TargetCalculatorTest::test_get_refmesh_Jacobian_2D()
01249 {
01250     MsqError err;
01251 
01252     const Vector3D tri_coords[] = { Vector3D( 0, 0, 0 ), Vector3D( 3, 0, 1 ), Vector3D( 1, 0, 3 ) };
01253     DummyRefMesh ref_mesh( tri_coords, 3 );
01254 
01255     const Vector3D right_coords[] = { Vector3D( 0, 0, 0 ), Vector3D( 5, 0, 0 ), Vector3D( 5, 1, 0 ) };
01256     size_t conn[]                 = { 0, 1, 2 };
01257     bool fixed[3];
01258     std::fill( fixed, fixed + sizeof( fixed ) / sizeof( fixed[0] ), false );
01259     PatchData pd;
01260     pd.fill( 3, right_coords[0].to_array(), 1, TRIANGLE, conn, fixed, err );
01261     ASSERT_NO_ERROR( err );
01262 
01263     Sample pts[]      = { Sample( 0, 0 ), Sample( 1, 1 ), Sample( 2, 0 ) };
01264     const int num_pts = sizeof( pts ) / sizeof( pts[0] );
01265     for( int i = 0; i < num_pts; ++i )
01266     {
01267         MsqMatrix< 3, 2 > J, W;
01268         TargetCalculator::jacobian_2D( pd, TRIANGLE, 3, pts[i], tri_coords, J, err );
01269         ASSERT_NO_ERROR( err );
01270 
01271         TargetCalculator::get_refmesh_Jacobian_2D( &ref_mesh, pd, 0, pts[i], W, err );
01272         ASSERT_NO_ERROR( err );
01273 
01274         ASSERT_MATRICES_EQUAL( J, W, EPS );
01275     }
01276 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines