MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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 }