MOAB: Mesh Oriented datABase  (version 5.2.1)
LinearMappingFunctionTest.cpp
Go to the documentation of this file.
00001 /* *****************************************************************
00002     MESQUITE -- The Mesh Quality Improvement Toolkit
00003 
00004     Copyright 2006 Lawrence Livermore National Laboratory.  Under
00005     the terms of Contract B545069 with the University of Wisconsin --
00006     Madison, Lawrence Livermore National Laboratory 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     (2006) kraftche@cae.wisc.edu
00024 
00025   ***************************************************************** */
00026 
00027 /** \file LinearMappingFunctionTest.cpp
00028  *  \brief
00029  *  \author Jason Kraftcheck
00030  */
00031 
00032 #include "Mesquite.hpp"
00033 #include "MappingFunction.hpp"
00034 #include "LinearHexahedron.hpp"
00035 #include "LinearQuadrilateral.hpp"
00036 #include "LinearTetrahedron.hpp"
00037 #include "LinearTriangle.hpp"
00038 #include "LinearPrism.hpp"
00039 #include "LinearPyramid.hpp"
00040 #include "TopologyInfo.hpp"
00041 #include "MsqError.hpp"
00042 #include "JacobianCalculator.hpp"
00043 #include "IdealElements.hpp"
00044 
00045 #include "UnitUtil.hpp"
00046 
00047 #include <vector>
00048 #include <algorithm>
00049 
00050 using namespace MBMesquite;
00051 using namespace std;
00052 
00053 class LinearMappingFunctionTest : public CppUnit::TestFixture
00054 {
00055   private:
00056     CPPUNIT_TEST_SUITE( LinearMappingFunctionTest );
00057 
00058     CPPUNIT_TEST( test_linear_hex_coeff_corners );
00059     CPPUNIT_TEST( test_linear_hex_coeff_edges );
00060     CPPUNIT_TEST( test_linear_hex_coeff_faces );
00061     CPPUNIT_TEST( test_linear_hex_coeff_center );
00062 
00063     CPPUNIT_TEST( test_linear_quad_coeff_corners );
00064     CPPUNIT_TEST( test_linear_quad_coeff_edges );
00065     //    CPPUNIT_TEST(test_linear_quad_coeff_faces);
00066     CPPUNIT_TEST( test_linear_quad_coeff_center );
00067 
00068     CPPUNIT_TEST( test_linear_tet_coeff_corners );
00069     CPPUNIT_TEST( test_linear_tet_coeff_edges );
00070     CPPUNIT_TEST( test_linear_tet_coeff_faces );
00071     CPPUNIT_TEST( test_linear_tet_coeff_center );
00072 
00073     CPPUNIT_TEST( test_linear_tri_coeff_corners );
00074     CPPUNIT_TEST( test_linear_tri_coeff_edges );
00075     CPPUNIT_TEST( test_linear_tri_coeff_faces );
00076     CPPUNIT_TEST( test_linear_tri_coeff_center );
00077 
00078     CPPUNIT_TEST( test_linear_prism_coeff_corners );
00079     CPPUNIT_TEST( test_linear_prism_coeff_edges );
00080     CPPUNIT_TEST( test_linear_prism_coeff_faces );
00081     CPPUNIT_TEST( test_linear_prism_coeff_center );
00082 
00083     CPPUNIT_TEST( test_linear_pyr_coeff_corners );
00084     CPPUNIT_TEST( test_linear_pyr_coeff_edges );
00085     CPPUNIT_TEST( test_linear_pyr_coeff_faces );
00086     CPPUNIT_TEST( test_linear_pyr_coeff_center );
00087 
00088     CPPUNIT_TEST( test_linear_hex_deriv_corners );
00089     CPPUNIT_TEST( test_linear_hex_deriv_edges );
00090     CPPUNIT_TEST( test_linear_hex_deriv_faces );
00091     CPPUNIT_TEST( test_linear_hex_deriv_center );
00092 
00093     CPPUNIT_TEST( test_linear_quad_deriv_corners );
00094     CPPUNIT_TEST( test_linear_quad_deriv_edges );
00095     //    CPPUNIT_TEST(test_linear_quad_deriv_faces);
00096     CPPUNIT_TEST( test_linear_quad_deriv_center );
00097 
00098     CPPUNIT_TEST( test_linear_tet_deriv_corners );
00099     CPPUNIT_TEST( test_linear_tet_deriv_edges );
00100     CPPUNIT_TEST( test_linear_tet_deriv_faces );
00101     CPPUNIT_TEST( test_linear_tet_deriv_center );
00102 
00103     CPPUNIT_TEST( test_linear_tri_deriv_corners );
00104     CPPUNIT_TEST( test_linear_tri_deriv_edges );
00105     //    CPPUNIT_TEST(test_linear_tri_deriv_faces);
00106     CPPUNIT_TEST( test_linear_tri_deriv_center );
00107 
00108     CPPUNIT_TEST( test_linear_prism_deriv_corners );
00109     CPPUNIT_TEST( test_linear_prism_deriv_edges );
00110     CPPUNIT_TEST( test_linear_prism_deriv_faces );
00111     CPPUNIT_TEST( test_linear_prism_deriv_center );
00112 
00113     CPPUNIT_TEST( test_linear_pyr_deriv_corners );
00114     CPPUNIT_TEST( test_linear_pyr_deriv_edges );
00115     CPPUNIT_TEST( test_linear_pyr_deriv_faces );
00116     CPPUNIT_TEST( test_linear_pyr_deriv_center );
00117 
00118     CPPUNIT_TEST( test_linear_hex_ideal );
00119     CPPUNIT_TEST( test_linear_quad_ideal );
00120     CPPUNIT_TEST( test_linear_tet_ideal );
00121     CPPUNIT_TEST( test_linear_tri_ideal );
00122     CPPUNIT_TEST( test_linear_prism_ideal );
00123 
00124     CPPUNIT_TEST_SUITE_END();
00125 
00126     LinearHexahedron hex;
00127     LinearQuadrilateral quad;
00128     LinearTetrahedron tet;
00129     LinearTriangle tri;
00130     LinearPrism prism;
00131     LinearPyramid pyr;
00132 
00133     static void hex_coeff( double xi[3], double coeff[8] );
00134     static void tet_coeff( double xi[3], double coeff[4] );
00135     static void quad_coeff( double xi[2], double coeff[4] );
00136     static void tri_coeff( double xi[2], double coeff[3] );
00137     static void prism_coeff( double xi[3], double coeff[6] );
00138     static void pyr_coeff( double xi[3], double coeff[5] );
00139 
00140     static void hex_deriv( double xi[3], double coeff_deriv[24] );
00141     static void tet_deriv( double xi[3], double coeff_deriv[12] );
00142     static void quad_deriv( double xi[2], double coeff_deriv[8] );
00143     static void tri_deriv( double xi[2], double coeff_deriv[6] );
00144     static void prism_deriv( double xi[3], double coeff_deriv[18] );
00145     static void pyr_deriv( double xi[3], double coeff_deriv[15] );
00146 
00147     typedef void ( *map_func )( double*, double* );
00148 
00149     void do_coeff_test( MappingFunction& mf, unsigned subdim, map_func mf2, unsigned count, double* xi );
00150     void do_deriv_test( MappingFunction2D& mf, unsigned subdim, map_func mf2, unsigned count, double* xi );
00151     void do_deriv_test( MappingFunction3D& mf, unsigned subdim, map_func mf2, unsigned count, double* xi );
00152     void do_ideal_test( MappingFunction2D& mf );
00153     void do_ideal_test( MappingFunction3D& mf );
00154 
00155     void test_coeff_fail( MappingFunction& mf, unsigned subdim );
00156     void test_deriv_fail( MappingFunction2D& mf, unsigned subdim );
00157     void test_deriv_fail( MappingFunction3D& mf, unsigned subdim );
00158 
00159     void xi_at_corners( EntityTopology type, double* xi, const int* corners );
00160     void xi_at_edges( EntityTopology type, double* xi, const int* corners );
00161     void xi_at_faces( EntityTopology type, double* xi, const int* corners );
00162 
00163   public:
00164     void setUp();
00165     void tearDown();
00166 
00167     void test_linear_hex_coeff_corners();
00168     void test_linear_hex_coeff_edges();
00169     void test_linear_hex_coeff_faces();
00170     void test_linear_hex_coeff_center();
00171 
00172     void test_linear_quad_coeff_corners();
00173     void test_linear_quad_coeff_edges();
00174     //    void test_linear_quad_coeff_faces();
00175     void test_linear_quad_coeff_center();
00176 
00177     void test_linear_tet_coeff_corners();
00178     void test_linear_tet_coeff_edges();
00179     void test_linear_tet_coeff_faces();
00180     void test_linear_tet_coeff_center();
00181 
00182     void test_linear_tri_coeff_corners();
00183     void test_linear_tri_coeff_edges();
00184     void test_linear_tri_coeff_faces();
00185     void test_linear_tri_coeff_center();
00186 
00187     void test_linear_prism_coeff_corners();
00188     void test_linear_prism_coeff_edges();
00189     void test_linear_prism_coeff_faces();
00190     void test_linear_prism_coeff_center();
00191 
00192     void test_linear_pyr_coeff_corners();
00193     void test_linear_pyr_coeff_edges();
00194     void test_linear_pyr_coeff_faces();
00195     void test_linear_pyr_coeff_center();
00196 
00197     void test_linear_hex_deriv_corners();
00198     void test_linear_hex_deriv_edges();
00199     void test_linear_hex_deriv_faces();
00200     void test_linear_hex_deriv_center();
00201 
00202     void test_linear_quad_deriv_corners();
00203     void test_linear_quad_deriv_edges();
00204     //    void test_linear_quad_deriv_faces();
00205     void test_linear_quad_deriv_center();
00206 
00207     void test_linear_tet_deriv_corners();
00208     void test_linear_tet_deriv_edges();
00209     void test_linear_tet_deriv_faces();
00210     void test_linear_tet_deriv_center();
00211 
00212     void test_linear_tri_deriv_corners();
00213     void test_linear_tri_deriv_edges();
00214     //    void test_linear_tri_deriv_faces();
00215     void test_linear_tri_deriv_center();
00216 
00217     void test_linear_prism_deriv_corners();
00218     void test_linear_prism_deriv_edges();
00219     void test_linear_prism_deriv_faces();
00220     void test_linear_prism_deriv_center();
00221 
00222     void test_linear_pyr_deriv_corners();
00223     void test_linear_pyr_deriv_edges();
00224     void test_linear_pyr_deriv_faces();
00225     void test_linear_pyr_deriv_center();
00226 
00227     void test_linear_hex_ideal()
00228     {
00229         do_ideal_test( hex );
00230     }
00231     void test_linear_quad_ideal()
00232     {
00233         do_ideal_test( quad );
00234     }
00235     void test_linear_tet_ideal()
00236     {
00237         do_ideal_test( tet );
00238     }
00239     void test_linear_tri_ideal()
00240     {
00241         do_ideal_test( tri );
00242     }
00243     void test_linear_prism_ideal()
00244     {
00245         do_ideal_test( prism );
00246     }
00247 };
00248 
00249 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( LinearMappingFunctionTest, "LinearMappingFunctionTest" );
00250 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( LinearMappingFunctionTest, "Unit" );
00251 
00252 void LinearMappingFunctionTest::setUp() {}
00253 void LinearMappingFunctionTest::tearDown() {}
00254 
00255 /*******************************************************************************
00256  *                             Xi values at element corners
00257  *******************************************************************************/
00258 
00259 const int HexCorners[8][3] = { { 0, 0, 0 }, { 1, 0, 0 }, { 1, 1, 0 }, { 0, 1, 0 },
00260                                { 0, 0, 1 }, { 1, 0, 1 }, { 1, 1, 1 }, { 0, 1, 1 } };
00261 
00262 const int QuadCorners[4][2] = { { 0, 0 }, { 1, 0 }, { 1, 1 }, { 0, 1 } };
00263 
00264 const int TetCorners[12] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1 };
00265 
00266 const int TriCorners[6] = { 0, 0, 1, 0, 0, 1 };
00267 
00268 const int PrismCorners[18] = { 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 1, 1, 0, 1, 0, 1 };
00269 
00270 /*******************************************************************************
00271  * Functions to calculate element xi at different locations, given the xi
00272  * at the corners.
00273  *******************************************************************************/
00274 
00275 void LinearMappingFunctionTest::xi_at_corners( EntityTopology type, double* xi, const int* corners )
00276 {
00277     unsigned d = TopologyInfo::dimension( type );
00278     unsigned c = TopologyInfo::corners( type );
00279     for( unsigned i = 0; i < c * d; ++i )
00280         xi[i] = corners[i];
00281 }
00282 
00283 void LinearMappingFunctionTest::xi_at_edges( EntityTopology type, double* xi, const int* corners )
00284 {
00285     MsqError err;
00286     unsigned d = TopologyInfo::dimension( type );
00287     unsigned e = TopologyInfo::edges( type );
00288     for( unsigned i = 0; i < e; ++i )
00289     {
00290         const unsigned* vtx = TopologyInfo::edge_vertices( type, i, err );
00291         CPPUNIT_ASSERT( !err );
00292         for( unsigned j = 0; j < d; ++j )
00293             xi[d * i + j] = ( corners[d * vtx[0] + j] + corners[d * vtx[1] + j] ) / 2.0;
00294     }
00295 }
00296 
00297 void LinearMappingFunctionTest::xi_at_faces( EntityTopology type, double* xi, const int* corners )
00298 {
00299     MsqError err;
00300     unsigned d = TopologyInfo::dimension( type );
00301     unsigned f = TopologyInfo::faces( type );
00302     for( unsigned i = 0; i < f; ++i )
00303     {
00304         unsigned c;
00305         const unsigned* vtx = TopologyInfo::face_vertices( type, i, c, err );
00306         CPPUNIT_ASSERT( !err );
00307         for( unsigned j = 0; j < d; ++j )
00308         {
00309             int sum = 0;
00310             for( unsigned k = 0; k < c; ++k )
00311                 sum += corners[d * vtx[k] + j];
00312             xi[d * i + j] = (double)sum / c;
00313         }
00314     }
00315 }
00316 
00317 /*******************************************************************************
00318  *                 Test functions for mapping function coefficients
00319  *******************************************************************************/
00320 
00321 void LinearMappingFunctionTest::test_linear_hex_coeff_corners()
00322 {
00323     double xi[24];
00324     xi_at_corners( HEXAHEDRON, xi, &HexCorners[0][0] );
00325     do_coeff_test( hex, 0, hex_coeff, 8, xi );
00326 }
00327 
00328 void LinearMappingFunctionTest::test_linear_hex_coeff_edges()
00329 {
00330     double xi[36];
00331     xi_at_edges( HEXAHEDRON, xi, &HexCorners[0][0] );
00332     do_coeff_test( hex, 1, hex_coeff, 12, xi );
00333 }
00334 
00335 void LinearMappingFunctionTest::test_linear_hex_coeff_faces()
00336 {
00337     double xi[18];
00338     xi_at_faces( HEXAHEDRON, xi, &HexCorners[0][0] );
00339     do_coeff_test( hex, 2, hex_coeff, 6, xi );
00340 }
00341 
00342 void LinearMappingFunctionTest::test_linear_hex_coeff_center()
00343 {
00344     double xi[3] = { 0.5, 0.5, 0.5 };
00345     do_coeff_test( hex, 3, hex_coeff, 1, xi );
00346 }
00347 
00348 void LinearMappingFunctionTest::test_linear_quad_coeff_corners()
00349 {
00350     double xi[8];
00351     xi_at_corners( QUADRILATERAL, xi, &QuadCorners[0][0] );
00352     do_coeff_test( quad, 0, quad_coeff, 4, xi );
00353 }
00354 
00355 void LinearMappingFunctionTest::test_linear_quad_coeff_edges()
00356 {
00357     double xi[8];
00358     xi_at_edges( QUADRILATERAL, xi, &QuadCorners[0][0] );
00359     do_coeff_test( quad, 1, quad_coeff, 4, xi );
00360 }
00361 
00362 // void LinearMappingFunctionTest::test_linear_quad_coeff_faces()
00363 //{
00364 //  test_coeff_fail( quad, 3 );
00365 //}
00366 
00367 void LinearMappingFunctionTest::test_linear_quad_coeff_center()
00368 {
00369     double xi[2] = { 0.5, 0.5 };
00370     do_coeff_test( quad, 2, quad_coeff, 1, xi );
00371 }
00372 
00373 void LinearMappingFunctionTest::test_linear_tet_coeff_corners()
00374 {
00375     double xi[12];
00376     xi_at_corners( TETRAHEDRON, xi, TetCorners );
00377     do_coeff_test( tet, 0, tet_coeff, 4, xi );
00378 }
00379 
00380 void LinearMappingFunctionTest::test_linear_tet_coeff_edges()
00381 {
00382     double xi[18];
00383     xi_at_edges( TETRAHEDRON, xi, TetCorners );
00384     do_coeff_test( tet, 1, tet_coeff, 6, xi );
00385 }
00386 
00387 void LinearMappingFunctionTest::test_linear_tet_coeff_faces()
00388 {
00389     double xi[12];
00390     xi_at_faces( TETRAHEDRON, xi, TetCorners );
00391     do_coeff_test( tet, 2, tet_coeff, 4, xi );
00392 }
00393 
00394 void LinearMappingFunctionTest::test_linear_tet_coeff_center()
00395 {
00396     double xi[3] = { 0.25, 0.25, 0.25 };
00397     do_coeff_test( tet, 3, tet_coeff, 1, xi );
00398 }
00399 
00400 void LinearMappingFunctionTest::test_linear_tri_coeff_corners()
00401 {
00402     double xi[12];
00403     xi_at_corners( TRIANGLE, xi, TriCorners );
00404     do_coeff_test( tri, 0, tri_coeff, 3, xi );
00405 }
00406 
00407 void LinearMappingFunctionTest::test_linear_tri_coeff_edges()
00408 {
00409     double xi[18];
00410     xi_at_edges( TRIANGLE, xi, TriCorners );
00411     do_coeff_test( tri, 1, tri_coeff, 3, xi );
00412 }
00413 
00414 void LinearMappingFunctionTest::test_linear_tri_coeff_faces()
00415 {
00416     test_coeff_fail( tri, 3 );
00417 }
00418 
00419 void LinearMappingFunctionTest::test_linear_tri_coeff_center()
00420 {
00421     double xi[2] = { 1. / 3, 1. / 3 };
00422     do_coeff_test( tri, 2, tri_coeff, 1, xi );
00423 }
00424 
00425 void LinearMappingFunctionTest::test_linear_prism_coeff_corners()
00426 {
00427     double xi[18];
00428     xi_at_corners( PRISM, xi, PrismCorners );
00429     do_coeff_test( prism, 0, prism_coeff, 6, xi );
00430 }
00431 
00432 void LinearMappingFunctionTest::test_linear_prism_coeff_edges()
00433 {
00434     double xi[27];
00435     xi_at_edges( PRISM, xi, PrismCorners );
00436     do_coeff_test( prism, 1, prism_coeff, 9, xi );
00437 }
00438 
00439 void LinearMappingFunctionTest::test_linear_prism_coeff_faces()
00440 {
00441     double xi[15];
00442     xi_at_faces( PRISM, xi, PrismCorners );
00443     do_coeff_test( prism, 2, prism_coeff, 5, xi );
00444 }
00445 
00446 void LinearMappingFunctionTest::test_linear_prism_coeff_center()
00447 {
00448     double xi[3] = { 0.5, 1. / 3, 1. / 3 };
00449     do_coeff_test( prism, 3, prism_coeff, 1, xi );
00450 }
00451 
00452 void LinearMappingFunctionTest::test_linear_pyr_coeff_corners()
00453 {
00454     double xi[15] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 1 };
00455     do_coeff_test( pyr, 0, pyr_coeff, 5, xi );
00456 }
00457 
00458 void LinearMappingFunctionTest::test_linear_pyr_coeff_edges()
00459 {
00460     double xi[24] = { 0.5, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.5, 0.0,
00461                       0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 1.0, 0.5, 0.0, 1.0, 0.5 };
00462     do_coeff_test( pyr, 1, pyr_coeff, 8, xi );
00463 }
00464 
00465 void LinearMappingFunctionTest::test_linear_pyr_coeff_faces()
00466 {
00467     double xi[15] = { 0.5, 0.0, 0.5, 1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.0 };
00468     do_coeff_test( pyr, 2, pyr_coeff, 5, xi );
00469 }
00470 
00471 void LinearMappingFunctionTest::test_linear_pyr_coeff_center()
00472 {
00473     double xi[3] = { 0.5, 0.5, 0.5 };
00474     do_coeff_test( pyr, 3, pyr_coeff, 1, xi );
00475 }
00476 
00477 /*******************************************************************************
00478  *                 Test functions for mapping function derivatives
00479  *******************************************************************************/
00480 
00481 void LinearMappingFunctionTest::test_linear_hex_deriv_corners()
00482 {
00483     double xi[24];
00484     xi_at_corners( HEXAHEDRON, xi, &HexCorners[0][0] );
00485     do_deriv_test( hex, 0, hex_deriv, 8, xi );
00486 }
00487 
00488 void LinearMappingFunctionTest::test_linear_hex_deriv_edges()
00489 {
00490     double xi[36];
00491     xi_at_edges( HEXAHEDRON, xi, &HexCorners[0][0] );
00492     do_deriv_test( hex, 1, hex_deriv, 12, xi );
00493 }
00494 
00495 void LinearMappingFunctionTest::test_linear_hex_deriv_faces()
00496 {
00497     double xi[18];
00498     xi_at_faces( HEXAHEDRON, xi, &HexCorners[0][0] );
00499     do_deriv_test( hex, 2, hex_deriv, 6, xi );
00500 }
00501 
00502 void LinearMappingFunctionTest::test_linear_hex_deriv_center()
00503 {
00504     double xi[3] = { 0.5, 0.5, 0.5 };
00505     do_deriv_test( hex, 3, hex_deriv, 1, xi );
00506 }
00507 
00508 void LinearMappingFunctionTest::test_linear_quad_deriv_corners()
00509 {
00510     double xi[8];
00511     xi_at_corners( QUADRILATERAL, xi, &QuadCorners[0][0] );
00512     do_deriv_test( quad, 0, quad_deriv, 4, xi );
00513 }
00514 
00515 void LinearMappingFunctionTest::test_linear_quad_deriv_edges()
00516 {
00517     double xi[8];
00518     xi_at_edges( QUADRILATERAL, xi, &QuadCorners[0][0] );
00519     do_deriv_test( quad, 1, quad_deriv, 4, xi );
00520 }
00521 
00522 // void LinearMappingFunctionTest::test_linear_quad_deriv_faces()
00523 //{
00524 //  test_deriv_fail( quad, 3 );
00525 //}
00526 
00527 void LinearMappingFunctionTest::test_linear_quad_deriv_center()
00528 {
00529     double xi[2] = { 0.5, 0.5 };
00530     do_deriv_test( quad, 2, quad_deriv, 1, xi );
00531 }
00532 
00533 void LinearMappingFunctionTest::test_linear_tet_deriv_corners()
00534 {
00535     double xi[12];
00536     xi_at_corners( TETRAHEDRON, xi, TetCorners );
00537     do_deriv_test( tet, 0, tet_deriv, 4, xi );
00538 }
00539 
00540 void LinearMappingFunctionTest::test_linear_tet_deriv_edges()
00541 {
00542     double xi[18];
00543     xi_at_edges( TETRAHEDRON, xi, TetCorners );
00544     do_deriv_test( tet, 1, tet_deriv, 6, xi );
00545 }
00546 
00547 void LinearMappingFunctionTest::test_linear_tet_deriv_faces()
00548 {
00549     double xi[12];
00550     xi_at_faces( TETRAHEDRON, xi, TetCorners );
00551     do_deriv_test( tet, 2, tet_deriv, 4, xi );
00552 }
00553 
00554 void LinearMappingFunctionTest::test_linear_tet_deriv_center()
00555 {
00556     double xi[3] = { 0.25, 0.25, 0.25 };
00557     do_deriv_test( tet, 3, tet_deriv, 1, xi );
00558 }
00559 
00560 void LinearMappingFunctionTest::test_linear_tri_deriv_corners()
00561 {
00562     double xi[12];
00563     xi_at_corners( TRIANGLE, xi, TriCorners );
00564     do_deriv_test( tri, 0, tri_deriv, 3, xi );
00565 }
00566 
00567 void LinearMappingFunctionTest::test_linear_tri_deriv_edges()
00568 {
00569     double xi[18];
00570     xi_at_edges( TRIANGLE, xi, TriCorners );
00571     do_deriv_test( tri, 1, tri_deriv, 3, xi );
00572 }
00573 
00574 // void LinearMappingFunctionTest::test_linear_tri_deriv_faces()
00575 //{
00576 //  test_deriv_fail( tri, 3 );
00577 //}
00578 
00579 void LinearMappingFunctionTest::test_linear_tri_deriv_center()
00580 {
00581     double xi[2] = { 1. / 3, 1. / 3 };
00582     do_deriv_test( tri, 2, tri_deriv, 1, xi );
00583 }
00584 
00585 void LinearMappingFunctionTest::test_linear_prism_deriv_corners()
00586 {
00587     double xi[18];
00588     xi_at_corners( PRISM, xi, PrismCorners );
00589     do_deriv_test( prism, 0, prism_deriv, 6, xi );
00590 }
00591 
00592 void LinearMappingFunctionTest::test_linear_prism_deriv_edges()
00593 {
00594     double xi[27];
00595     xi_at_edges( PRISM, xi, PrismCorners );
00596     do_deriv_test( prism, 1, prism_deriv, 9, xi );
00597 }
00598 
00599 void LinearMappingFunctionTest::test_linear_prism_deriv_faces()
00600 {
00601     double xi[15];
00602     xi_at_faces( PRISM, xi, PrismCorners );
00603     do_deriv_test( prism, 2, prism_deriv, 5, xi );
00604 }
00605 
00606 void LinearMappingFunctionTest::test_linear_prism_deriv_center()
00607 {
00608     double xi[3] = { 0.5, 1. / 3, 1. / 3 };
00609     do_deriv_test( prism, 3, prism_deriv, 1, xi );
00610 }
00611 
00612 void LinearMappingFunctionTest::test_linear_pyr_deriv_corners()
00613 {
00614     double xi[15] = { 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0.5, 0.5, 1 };
00615     do_deriv_test( pyr, 0, pyr_deriv, 5, xi );
00616 }
00617 
00618 void LinearMappingFunctionTest::test_linear_pyr_deriv_edges()
00619 {
00620     double xi[24] = { 0.5, 0.0, 0.0, 1.0, 0.5, 0.0, 0.5, 1.0, 0.0, 0.0, 0.5, 0.0,
00621                       0.0, 0.0, 0.5, 1.0, 0.0, 0.5, 1.0, 1.0, 0.5, 0.0, 1.0, 0.5 };
00622     do_deriv_test( pyr, 1, pyr_deriv, 8, xi );
00623 }
00624 
00625 void LinearMappingFunctionTest::test_linear_pyr_deriv_faces()
00626 {
00627     double xi[15] = { 0.5, 0.0, 0.5, 1.0, 0.5, 0.5, 0.5, 1.0, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.0 };
00628     do_deriv_test( pyr, 2, pyr_deriv, 5, xi );
00629 }
00630 
00631 void LinearMappingFunctionTest::test_linear_pyr_deriv_center()
00632 {
00633     double xi[3] = { 0.5, 0.5, 0.5 };
00634     do_deriv_test( pyr, 3, pyr_deriv, 1, xi );
00635 }
00636 
00637 /*******************************************************************************
00638  *        Mapping function implementations to compare with
00639  *******************************************************************************/
00640 
00641 void LinearMappingFunctionTest::hex_coeff( double xi[3], double coeff[8] )
00642 {
00643     coeff[0] = ( 1 - xi[0] ) * ( 1 - xi[1] ) * ( 1 - xi[2] );
00644     coeff[1] = xi[0] * ( 1 - xi[1] ) * ( 1 - xi[2] );
00645     coeff[2] = xi[0] * xi[1] * ( 1 - xi[2] );
00646     coeff[3] = ( 1 - xi[0] ) * xi[1] * ( 1 - xi[2] );
00647     coeff[4] = ( 1 - xi[0] ) * ( 1 - xi[1] ) * xi[2];
00648     coeff[5] = xi[0] * ( 1 - xi[1] ) * xi[2];
00649     coeff[6] = xi[0] * xi[1] * xi[2];
00650     coeff[7] = ( 1 - xi[0] ) * xi[1] * xi[2];
00651 }
00652 
00653 void LinearMappingFunctionTest::tet_coeff( double xi[3], double coeff[4] )
00654 {
00655     coeff[0] = 1 - xi[0] - xi[1] - xi[2];
00656     coeff[1] = xi[0];
00657     coeff[2] = xi[1];
00658     coeff[3] = xi[2];
00659 }
00660 
00661 void LinearMappingFunctionTest::quad_coeff( double xi[2], double coeff[4] )
00662 {
00663     coeff[0] = ( 1 - xi[0] ) * ( 1 - xi[1] );
00664     coeff[1] = xi[0] * ( 1 - xi[1] );
00665     coeff[2] = xi[0] * xi[1];
00666     coeff[3] = ( 1 - xi[0] ) * xi[1];
00667 }
00668 
00669 void LinearMappingFunctionTest::tri_coeff( double xi[2], double coeff[3] )
00670 {
00671     coeff[0] = 1 - xi[0] - xi[1];
00672     coeff[1] = xi[0];
00673     coeff[2] = xi[1];
00674 }
00675 
00676 void LinearMappingFunctionTest::prism_coeff( double xi[3], double coeff[6] )
00677 {
00678     coeff[0] = ( 1 - xi[0] ) * ( 1 - xi[1] - xi[2] );
00679     coeff[1] = ( 1 - xi[0] ) * xi[1];
00680     coeff[2] = ( 1 - xi[0] ) * xi[2];
00681     coeff[3] = xi[0] * ( 1 - xi[1] - xi[2] );
00682     coeff[4] = xi[0] * xi[1];
00683     coeff[5] = xi[0] * xi[2];
00684 }
00685 
00686 void LinearMappingFunctionTest::pyr_coeff( double xi[3], double coeff[5] )
00687 {
00688     coeff[0] = ( 1 - xi[0] ) * ( 1 - xi[1] ) * ( 1 - xi[2] );
00689     coeff[1] = xi[0] * ( 1 - xi[1] ) * ( 1 - xi[2] );
00690     coeff[2] = xi[0] * xi[1] * ( 1 - xi[2] );
00691     coeff[3] = ( 1 - xi[0] ) * xi[1] * ( 1 - xi[2] );
00692     coeff[4] = xi[2];
00693 }
00694 
00695 /*******************************************************************************
00696  *        Mapping function derivatives to compare with
00697  *******************************************************************************/
00698 
00699 void LinearMappingFunctionTest::hex_deriv( double xi[3], double coeff[24] )
00700 {
00701     coeff[3 * 0 + 0] = -( 1 - xi[1] ) * ( 1 - xi[2] );
00702     coeff[3 * 0 + 1] = -( 1 - xi[0] ) * ( 1 - xi[2] );
00703     coeff[3 * 0 + 2] = -( 1 - xi[0] ) * ( 1 - xi[1] );
00704 
00705     coeff[3 * 1 + 0] = ( 1 - xi[1] ) * ( 1 - xi[2] );
00706     coeff[3 * 1 + 1] = -xi[0] * ( 1 - xi[2] );
00707     coeff[3 * 1 + 2] = -xi[0] * ( 1 - xi[1] );
00708 
00709     coeff[3 * 2 + 0] = xi[1] * ( 1 - xi[2] );
00710     coeff[3 * 2 + 1] = xi[0] * ( 1 - xi[2] );
00711     coeff[3 * 2 + 2] = -xi[0] * xi[1];
00712 
00713     coeff[3 * 3 + 0] = -xi[1] * ( 1 - xi[2] );
00714     coeff[3 * 3 + 1] = ( 1 - xi[0] ) * ( 1 - xi[2] );
00715     coeff[3 * 3 + 2] = -( 1 - xi[0] ) * xi[1];
00716 
00717     coeff[3 * 4 + 0] = -( 1 - xi[1] ) * xi[2];
00718     coeff[3 * 4 + 1] = -( 1 - xi[0] ) * xi[2];
00719     coeff[3 * 4 + 2] = ( 1 - xi[0] ) * ( 1 - xi[1] );
00720 
00721     coeff[3 * 5 + 0] = ( 1 - xi[1] ) * xi[2];
00722     coeff[3 * 5 + 1] = -xi[0] * xi[2];
00723     coeff[3 * 5 + 2] = xi[0] * ( 1 - xi[1] );
00724 
00725     coeff[3 * 6 + 0] = xi[1] * xi[2];
00726     coeff[3 * 6 + 1] = xi[0] * xi[2];
00727     coeff[3 * 6 + 2] = xi[0] * xi[1];
00728 
00729     coeff[3 * 7 + 0] = -xi[1] * xi[2];
00730     coeff[3 * 7 + 1] = ( 1 - xi[0] ) * xi[2];
00731     coeff[3 * 7 + 2] = ( 1 - xi[0] ) * xi[1];
00732 }
00733 
00734 void LinearMappingFunctionTest::tet_deriv( double*, double coeff[12] )
00735 {
00736     static const double derivs[] = { -1, -1, -1, 1, 0, 0, 0, 1, 0, 0, 0, 1 };
00737     memcpy( coeff, derivs, sizeof( derivs ) );
00738 }
00739 
00740 void LinearMappingFunctionTest::quad_deriv( double xi[2], double coeff[8] )
00741 {
00742     coeff[2 * 0 + 0] = xi[1] - 1;
00743     coeff[2 * 0 + 1] = xi[0] - 1;
00744 
00745     coeff[2 * 1 + 0] = 1 - xi[1];
00746     coeff[2 * 1 + 1] = -xi[0];
00747 
00748     coeff[2 * 2 + 0] = xi[1];
00749     coeff[2 * 2 + 1] = xi[0];
00750 
00751     coeff[2 * 3 + 0] = -xi[1];
00752     coeff[2 * 3 + 1] = 1 - xi[0];
00753 }
00754 
00755 void LinearMappingFunctionTest::tri_deriv( double*, double coeff[6] )
00756 {
00757     static const double derivs[] = { -1, -1, 1, 0, 0, 1 };
00758     memcpy( coeff, derivs, sizeof( derivs ) );
00759 }
00760 
00761 void LinearMappingFunctionTest::prism_deriv( double xi[3], double coeff[18] )
00762 {
00763     coeff[0] = xi[1] + xi[2] - 1.0;
00764     ;
00765     coeff[1] = xi[0] - 1.0;
00766     coeff[2] = xi[0] - 1.0;
00767 
00768     coeff[3] = -xi[1];
00769     coeff[4] = 1.0 - xi[0];
00770     coeff[5] = 0.0;
00771 
00772     coeff[6] = -xi[2];
00773     coeff[7] = 0.0;
00774     coeff[8] = 1.0 - xi[0];
00775 
00776     coeff[9]  = 1.0 - xi[1] - xi[2];
00777     coeff[10] = -xi[0];
00778     coeff[11] = -xi[0];
00779 
00780     coeff[12] = xi[1];
00781     coeff[13] = xi[0];
00782     coeff[14] = 0.0;
00783 
00784     coeff[15] = xi[2];
00785     coeff[16] = 0.0;
00786     coeff[17] = xi[0];
00787 }
00788 
00789 void LinearMappingFunctionTest::pyr_deriv( double xi[3], double coeff[15] )
00790 {
00791     coeff[3 * 0 + 0] = -( 1 - xi[1] ) * ( 1 - xi[2] );
00792     coeff[3 * 0 + 1] = -( 1 - xi[0] ) * ( 1 - xi[2] );
00793     coeff[3 * 0 + 2] = -( 1 - xi[0] ) * ( 1 - xi[1] );
00794 
00795     coeff[3 * 1 + 0] = ( 1 - xi[1] ) * ( 1 - xi[2] );
00796     coeff[3 * 1 + 1] = -xi[0] * ( 1 - xi[2] );
00797     coeff[3 * 1 + 2] = -xi[0] * ( 1 - xi[1] );
00798 
00799     coeff[3 * 2 + 0] = xi[1] * ( 1 - xi[2] );
00800     coeff[3 * 2 + 1] = xi[0] * ( 1 - xi[2] );
00801     coeff[3 * 2 + 2] = -xi[0] * xi[1];
00802 
00803     coeff[3 * 3 + 0] = -xi[1] * ( 1 - xi[2] );
00804     coeff[3 * 3 + 1] = ( 1 - xi[0] ) * ( 1 - xi[2] );
00805     coeff[3 * 3 + 2] = -xi[1] * ( 1 - xi[0] );
00806 
00807     coeff[3 * 4 + 0] = 0.0;
00808     coeff[3 * 4 + 1] = 0.0;
00809     coeff[3 * 4 + 2] = 1.0;
00810 }
00811 
00812 /*******************************************************************************
00813  *         Some utlity functions for creating CppUnit error messages
00814  *******************************************************************************/
00815 
00816 static string itostr( int i )
00817 {
00818     char buffer[32];
00819     sprintf( buffer, "%d", i );
00820     return buffer;
00821 }
00822 
00823 static string dtostr( double i )
00824 {
00825     char buffer[32];
00826     sprintf( buffer, "%g", i );
00827     return buffer;
00828 }
00829 
00830 /*******************************************************************************
00831  *         Actual test imlplementation (common code for many tests)
00832  *******************************************************************************/
00833 
00834 void LinearMappingFunctionTest::do_coeff_test( MappingFunction& mf, unsigned subdim, map_func mf2, unsigned count,
00835                                                double* xi )
00836 {
00837     // make sure it fails if passed a nonlinear element
00838     MsqError err;
00839     double coeff[100];
00840     size_t indices[100];
00841     size_t num_coeff = 100;
00842     NodeSet tmp_set;
00843     tmp_set.set_mid_edge_node( 1 );
00844     mf.coefficients( Sample( 0, 1 ), tmp_set, coeff, indices, num_coeff, err );
00845     CPPUNIT_ASSERT( err );
00846     err.clear();
00847 
00848     // get number of vertices in element
00849     const unsigned n = TopologyInfo::corners( mf.element_topology() );
00850     const unsigned d = TopologyInfo::dimension( mf.element_topology() );
00851 
00852     // compare coefficients at each location
00853     vector< double > comp( n );
00854     for( unsigned i = 0; i < count; ++i )
00855     {
00856         num_coeff = 101;
00857         mf.coefficients( Sample( subdim, i ), NodeSet(), coeff, indices, num_coeff, err );
00858         CPPUNIT_ASSERT( !err );
00859 
00860         mf2( xi, &comp[0] );
00861         string xi_str;
00862         for( unsigned j = 0; j < d; ++j )
00863         {
00864             xi_str += !j ? "(" : ", ";
00865             xi_str += dtostr( xi[j] );
00866         }
00867         xi_str += ")";
00868         xi += d;
00869 
00870         for( unsigned j = 0; j < n; ++j )
00871         {
00872             double mf_val = 0.0;
00873             size_t idx    = std::find( indices, indices + num_coeff, j ) - indices;
00874             if( idx < num_coeff ) mf_val = coeff[idx];
00875 
00876             CppUnit::Message message( "Coefficients do not match." );
00877             message.addDetail( string( "Entity:             " ) + itostr( i ) );
00878             message.addDetail( string( "Coefficient number: " ) + itostr( j ) );
00879             message.addDetail( string( "Xi:             " ) + xi_str );
00880             message.addDetail( string( "Expected value: " ) + dtostr( comp[j] ) );
00881             message.addDetail( string( "Actual value:   " ) + dtostr( mf_val ) );
00882             ASSERT_MESSAGE( message, fabs( comp[j] - mf_val ) < DBL_EPSILON );
00883         }
00884     }
00885 }
00886 
00887 void LinearMappingFunctionTest::do_deriv_test( MappingFunction2D& mf, unsigned subdim, map_func mf2, unsigned count,
00888                                                double* xi )
00889 {
00890     // make sure it fails if passed a nonlinear element
00891     MsqError err;
00892     MsqVector< 2 > derivs[100];
00893     size_t verts[100], num_vtx = 37;
00894     NodeSet tmp_set;
00895     tmp_set.set_mid_edge_node( 1 );
00896     mf.derivatives( Sample( subdim, 0 ), tmp_set, verts, derivs, num_vtx, err );
00897     CPPUNIT_ASSERT( err );
00898     err.clear();
00899 
00900     // get number of vertices in element
00901     const unsigned n = TopologyInfo::corners( mf.element_topology() );
00902 
00903     // compare coefficients at each location
00904     vector< double > comp( 2 * n );
00905     for( unsigned i = 0; i < count; ++i )
00906     {
00907         num_vtx = 33;
00908         mf.derivatives( Sample( subdim, i ), NodeSet(), verts, derivs, num_vtx, err );
00909         CPPUNIT_ASSERT( !err );
00910         CPPUNIT_ASSERT( num_vtx > 0 );
00911         CPPUNIT_ASSERT( num_vtx <= n );
00912 
00913         mf2( xi, &comp[0] );
00914         string xi_str;
00915         for( unsigned j = 0; j < 2; ++j )
00916         {
00917             xi_str += !j ? "(" : ", ";
00918             xi_str += dtostr( xi[j] );
00919         }
00920         xi_str += ")";
00921         xi += 2;
00922 
00923         for( unsigned j = 0; j < num_vtx; ++j )
00924         {
00925             bool all_zero = true;
00926             for( unsigned k = 0; k < 2; ++k )
00927             {
00928                 CppUnit::Message message( "Coefficient derivatives do not match." );
00929                 message.addDetail( string( "Entity:             " ) + itostr( i ) );
00930                 message.addDetail( string( "Coefficient number: " ) + itostr( j ) );
00931                 message.addDetail( string( "Xi:             " ) + xi_str );
00932                 message.addDetail( string( "Axis:           " ) + itostr( k ) );
00933                 message.addDetail( string( "Expected value: " ) + dtostr( comp[2 * verts[j] + k] ) );
00934                 message.addDetail( string( "Actual value:   " ) + dtostr( derivs[j][k] ) );
00935                 ASSERT_MESSAGE( message, fabs( comp[2 * verts[j] + k] - derivs[j][k] ) < DBL_EPSILON );
00936                 if( fabs( derivs[j][k] ) > DBL_EPSILON ) all_zero = false;
00937             }
00938 
00939             // if vertex has all zero values, it shouldn't have been in the
00940             // vertex list at all, as the Jacobian will not depend on that vertex.
00941             CPPUNIT_ASSERT( !all_zero );
00942         }
00943 
00944         // If any vertex is not in the list, then its values must be zero.
00945         sort( verts, verts + num_vtx );
00946         for( unsigned j = 0; j < num_vtx; ++j )
00947         {
00948             if( !binary_search( verts, verts + num_vtx, j ) )
00949             {
00950                 for( unsigned k = 0; k < 2; ++k )
00951                 {
00952                     CppUnit::Message message( "Missing coefficient derivatives." );
00953                     message.addDetail( string( "Entity:              " ) + itostr( i ) );
00954                     message.addDetail( string( "Coefficient number:  " ) + itostr( j ) );
00955                     message.addDetail( string( "Axis:                " ) + itostr( k ) );
00956                     message.addDetail( string( "Expected derivative: " ) + dtostr( comp[2 * j + k] ) );
00957                     ASSERT_MESSAGE( message, fabs( comp[2 * j + k] ) < DBL_EPSILON );
00958                 }
00959             }
00960         }
00961     }
00962 }
00963 
00964 void LinearMappingFunctionTest::do_deriv_test( MappingFunction3D& mf, unsigned subdim, map_func mf2, unsigned count,
00965                                                double* xi )
00966 {
00967     // make sure it fails if passed a nonlinear element
00968     MsqError err;
00969     MsqVector< 3 > derivs[100];
00970     size_t verts[100], num_vtx = 37;
00971     NodeSet tmp_set;
00972     tmp_set.set_mid_edge_node( 1 );
00973     mf.derivatives( Sample( subdim, 0 ), tmp_set, verts, derivs, num_vtx, err );
00974     CPPUNIT_ASSERT( err );
00975     err.clear();
00976 
00977     // get number of vertices in element
00978     const unsigned n = TopologyInfo::corners( mf.element_topology() );
00979 
00980     // compare coefficients at each location
00981     vector< double > comp( 3 * n );
00982     for( unsigned i = 0; i < count; ++i )
00983     {
00984         num_vtx = 33;
00985         mf.derivatives( Sample( subdim, i ), NodeSet(), verts, derivs, num_vtx, err );
00986         CPPUNIT_ASSERT( !err );
00987         CPPUNIT_ASSERT( num_vtx > 0 );
00988         CPPUNIT_ASSERT( num_vtx <= n );
00989 
00990         mf2( xi, &comp[0] );
00991         string xi_str;
00992         for( unsigned j = 0; j < 3; ++j )
00993         {
00994             xi_str += !j ? "(" : ", ";
00995             xi_str += dtostr( xi[j] );
00996         }
00997         xi_str += ")";
00998         xi += 3;
00999 
01000         for( unsigned j = 0; j < num_vtx; ++j )
01001         {
01002             bool all_zero = true;
01003             for( unsigned k = 0; k < 3; ++k )
01004             {
01005                 CppUnit::Message message( "Coefficient derivatives do not match." );
01006                 message.addDetail( string( "Entity:             " ) + itostr( i ) );
01007                 message.addDetail( string( "Coefficient number: " ) + itostr( j ) );
01008                 message.addDetail( string( "Xi:             " ) + xi_str );
01009                 message.addDetail( string( "Axis:           " ) + itostr( k ) );
01010                 message.addDetail( string( "Expected value: " ) + dtostr( comp[3 * verts[j] + k] ) );
01011                 message.addDetail( string( "Actual value:   " ) + dtostr( derivs[j][k] ) );
01012                 ASSERT_MESSAGE( message, fabs( comp[3 * verts[j] + k] - derivs[j][k] ) < DBL_EPSILON );
01013                 if( fabs( derivs[j][k] ) > DBL_EPSILON ) all_zero = false;
01014             }
01015 
01016             // if vertex has all zero values, it shouldn't have been in the
01017             // vertex list at all, as the Jacobian will not depend on that vertex.
01018             CPPUNIT_ASSERT( !all_zero );
01019         }
01020 
01021         // If any vertex is not in the list, then its values must be zero.
01022         sort( verts, verts + num_vtx );
01023         for( unsigned j = 0; j < num_vtx; ++j )
01024         {
01025             if( !binary_search( verts, verts + num_vtx, j ) )
01026             {
01027                 for( unsigned k = 0; k < 3; ++k )
01028                 {
01029                     CppUnit::Message message( "Missing coefficient derivatives." );
01030                     message.addDetail( string( "Entity:              " ) + itostr( i ) );
01031                     message.addDetail( string( "Coefficient number:  " ) + itostr( j ) );
01032                     message.addDetail( string( "Axis:                " ) + itostr( k ) );
01033                     message.addDetail( string( "Expected derivative: " ) + dtostr( comp[3 * j + k] ) );
01034                     ASSERT_MESSAGE( message, fabs( comp[3 * j + k] ) < DBL_EPSILON );
01035                 }
01036             }
01037         }
01038     }
01039 }
01040 
01041 void LinearMappingFunctionTest::do_ideal_test( MappingFunction2D& mf )
01042 {
01043     MsqError err;
01044     MsqMatrix< 3, 2 > W_prime;
01045     mf.ideal( Sample( 2, 0 ), W_prime, err );
01046     ASSERT_NO_ERROR( err );
01047 
01048     // for this test that everything is in the xy-plane
01049     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 0 ), 1e-12 );
01050     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 1 ), 1e-12 );
01051     MsqMatrix< 2, 2 > W( W_prime.data() );
01052     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( W ), 1e-6 );
01053 
01054     const Vector3D* verts = unit_edge_element( mf.element_topology() );
01055     CPPUNIT_ASSERT( verts );
01056 
01057     JacobianCalculator jc;
01058     jc.get_Jacobian_2D( &mf, NodeSet(), Sample( 2, 0 ), verts, TopologyInfo::corners( mf.element_topology() ), W_prime,
01059                         err );
01060     ASSERT_NO_ERROR( err );
01061 
01062     // for this test that everything is in the xy-plane
01063     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 0 ), 1e-12 );
01064     CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.0, W_prime( 2, 1 ), 1e-12 );
01065     MsqMatrix< 2, 2 > W_exp( W_prime.data() );
01066     W_exp /= sqrt( det( W_exp ) );
01067 
01068     // Matrices should be a rotation of each other.
01069     // First, calculate tentative rotation matrix
01070     MsqMatrix< 2, 2 > R = inverse( W_exp ) * W;
01071     // next check that it is a rotation
01072     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );          // no scaling
01073     ASSERT_MATRICES_EQUAL( transpose( R ), inverse( R ), 1e-6 );  // orthogonal
01074 }
01075 
01076 void LinearMappingFunctionTest::do_ideal_test( MappingFunction3D& mf )
01077 {
01078     MsqError err;
01079     MsqMatrix< 3, 3 > W, I( 1.0 );
01080     mf.ideal( Sample( 3, 0 ), W, err );
01081     ASSERT_NO_ERROR( err );
01082     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( W ), 1e-6 );
01083 
01084     const Vector3D* verts = unit_edge_element( mf.element_topology() );
01085     CPPUNIT_ASSERT( verts );
01086 
01087     JacobianCalculator jc;
01088     MsqMatrix< 3, 3 > W_exp;
01089     jc.get_Jacobian_3D( &mf, NodeSet(), Sample( 3, 0 ), verts, TopologyInfo::corners( mf.element_topology() ), W_exp,
01090                         err );
01091     ASSERT_NO_ERROR( err );
01092     W_exp /= MBMesquite::cbrt( det( W_exp ) );
01093 
01094     // Matrices should be a rotation of each other.
01095     // First, calculate tentative rotation matrix
01096     MsqMatrix< 3, 3 > R = W * inverse( W_exp );
01097     // next check that it is a rotation
01098     CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, det( R ), 1e-6 );   // no scaling
01099     ASSERT_MATRICES_EQUAL( I, transpose( R ) * R, 1e-6 );  // orthogonal
01100 }
01101 
01102 void LinearMappingFunctionTest::test_coeff_fail( MappingFunction& mf, unsigned subdim )
01103 {
01104     // make sure it fails if called
01105     MsqError err;
01106     double coeff[100];
01107     size_t num_coeff, indices[100];
01108     mf.coefficients( Sample( subdim, 0 ), NodeSet(), coeff, indices, num_coeff, err );
01109     CPPUNIT_ASSERT( err );
01110     err.clear();
01111 }
01112 
01113 void LinearMappingFunctionTest::test_deriv_fail( MappingFunction2D& mf, unsigned subdim )
01114 {
01115     // make sure it fails if called
01116     MsqError err;
01117     MsqVector< 2 > coeff[100];
01118     size_t verts[100], num_coeff;
01119     mf.derivatives( Sample( subdim, 0 ), NodeSet(), verts, coeff, num_coeff, err );
01120     CPPUNIT_ASSERT( err );
01121     err.clear();
01122 }
01123 
01124 void LinearMappingFunctionTest::test_deriv_fail( MappingFunction3D& mf, unsigned subdim )
01125 {
01126     // make sure it fails if called
01127     MsqError err;
01128     MsqVector< 3 > coeff[100];
01129     size_t verts[100], num_coeff;
01130     mf.derivatives( Sample( subdim, 0 ), NodeSet(), verts, coeff, num_coeff, err );
01131     CPPUNIT_ASSERT( err );
01132     err.clear();
01133 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines