MOAB: Mesh Oriented datABase
(version 5.2.1)
|
00001 /* ***************************************************************** 00002 MESQUITE -- The Mesh Quality Improvement Toolkit 00003 00004 Copyright 2007 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 HigherOrderTest.cpp 00028 * \brief 00029 * \author Jason Kraftcheck 00030 */ 00031 00032 #define HAVE_HO_HEX 00033 #define TEST_HO_QUAD 00034 00035 #include "Mesquite.hpp" 00036 #include "MsqError.hpp" 00037 00038 #include "ArrayMesh.hpp" 00039 #include "DomainClassifier.hpp" 00040 #include "PlanarDomain.hpp" 00041 #include "MeshDomain1D.hpp" 00042 00043 #include "IdealShapeTarget.hpp" 00044 #include "TInverseMeanRatio.hpp" 00045 #include "TShapeSizeOrientNB1.hpp" 00046 #include "TShapeSizeNB3.hpp" 00047 #include "TQualityMetric.hpp" 00048 #include "PMeanPTemplate.hpp" 00049 #include "LPtoPTemplate.hpp" 00050 #include "SteepestDescent.hpp" 00051 #include "FeasibleNewton.hpp" 00052 #include "InstructionQueue.hpp" 00053 #include "TerminationCriterion.hpp" 00054 00055 //#include "TriLagrangeShape.hpp" 00056 //#include "TetLagrangeShape.hpp" 00057 //#include "QuadLagrangeShape.hpp" 00058 #ifdef HAVE_HO_HEX 00059 #include "HexLagrangeShape.hpp" 00060 #endif 00061 00062 #include "UnitUtil.hpp" 00063 00064 using namespace MBMesquite; 00065 00066 #include <iostream> 00067 using std::cerr; 00068 using std::cout; 00069 using std::endl; 00070 using std::vector; 00071 00072 const int MAX_ITERATIONS = 1000; 00073 const double QEL = 1.0; // edge length of ideal quad 00074 00075 class HigherOrderTest : public CppUnit::TestFixture 00076 { 00077 private: 00078 CPPUNIT_TEST_SUITE( HigherOrderTest ); 00079 CPPUNIT_TEST( test_tri_basic_ideal ); 00080 CPPUNIT_TEST( test_tri_basic_mid_spin ); 00081 CPPUNIT_TEST( test_tri_basic_mid_convex ); 00082 CPPUNIT_TEST( test_tri_basic_peak_up ); 00083 CPPUNIT_TEST( test_tri_basic_peak_down ); 00084 CPPUNIT_TEST( test_tri_basic_peak_over ); 00085 #ifdef TEST_HO_QUAD 00086 CPPUNIT_TEST( test_quad_basic_ideal ); 00087 CPPUNIT_TEST( test_quad_basic_mid_spin ); 00088 CPPUNIT_TEST( test_quad_basic_mid_convex ); 00089 CPPUNIT_TEST( test_quad_basic_left_down ); 00090 CPPUNIT_TEST( test_quad_basic_top_down ); 00091 CPPUNIT_TEST( test_quad_basic_right_up ); 00092 CPPUNIT_TEST( test_quad_basic_left_over ); 00093 #endif 00094 CPPUNIT_TEST( test_tet_basic_ideal ); 00095 CPPUNIT_TEST( test_tet_basic_mid_spin ); 00096 CPPUNIT_TEST( test_tet_basic_mid_convex ); 00097 CPPUNIT_TEST( test_tet_basic_apex_up ); 00098 CPPUNIT_TEST( test_tet_basic_apex_down ); 00099 CPPUNIT_TEST( test_tet_basic_apex_over ); 00100 #ifdef HAVE_HO_HEX 00101 CPPUNIT_TEST( test_hex_basic ); 00102 #endif 00103 // CPPUNIT_TEST (test_tri_open_domain); 00104 // CPPUNIT_TEST (test_tri_slac); 00105 CPPUNIT_TEST_SUITE_END(); 00106 00107 // TriLagrangeShape tri_shape; 00108 // TetLagrangeShape tet_shape; 00109 // QuadLagrangeShape quad_shape; 00110 #ifdef HAVE_HO_HEX 00111 HexLagrangeShape hex_shape; 00112 #endif 00113 InstructionQueue q; 00114 TShapeSizeNB3 tm; 00115 IdealShapeTarget tc; 00116 TQualityMetric metric; 00117 PMeanPTemplate func; 00118 SteepestDescent solver; 00119 TerminationCriterion crit, outer; 00120 00121 public: 00122 HigherOrderTest() : metric( &tc, &tm ), func( 1, &metric ), solver( &func ) 00123 { 00124 MsqError err; 00125 // q.set_mapping_function( &tri_shape ); 00126 // q.set_mapping_function( &tet_shape ); 00127 // q.set_mapping_function( &quad_shape ); 00128 #ifdef HAVE_HO_HEX 00129 q.set_mapping_function( &hex_shape ); 00130 #endif 00131 q.set_master_quality_improver( &solver, err ); 00132 00133 q.set_slaved_ho_node_mode( Settings::SLAVE_NONE ); 00134 00135 outer.add_iteration_limit( 1 ); 00136 crit.add_absolute_vertex_movement( 1e-6 ); 00137 crit.add_iteration_limit( MAX_ITERATIONS ); 00138 solver.set_outer_termination_criterion( &outer ); 00139 solver.set_inner_termination_criterion( &crit ); 00140 } 00141 00142 bool hit_iteration_limit() const 00143 { 00144 return crit.get_iteration_count() >= MAX_ITERATIONS; 00145 } 00146 00147 void test_tri_basic_ideal(); 00148 void test_tri_basic_mid_spin(); 00149 void test_tri_basic_mid_convex(); 00150 void test_tri_basic_peak_up(); 00151 void test_tri_basic_peak_down(); 00152 void test_tri_basic_peak_over(); 00153 void test_quad_basic_ideal(); 00154 void test_quad_basic_mid_spin(); 00155 void test_quad_basic_mid_convex(); 00156 void test_quad_basic_left_down(); 00157 void test_quad_basic_top_down(); 00158 void test_quad_basic_right_up(); 00159 void test_quad_basic_left_over(); 00160 void test_tet_basic_ideal(); 00161 void test_tet_basic_mid_spin(); 00162 void test_tet_basic_mid_convex(); 00163 void test_tet_basic_apex_up(); 00164 void test_tet_basic_apex_down(); 00165 void test_tet_basic_apex_over(); 00166 void test_hex_basic(); 00167 void test_tri_open_domain(); 00168 void test_tri_slac(); 00169 00170 void test_tri_open_domain( double& x1, double& x3, double& x4, double y2, double y5, double& y4 ); 00171 00172 void basic_tri_test( double& x2, double& y2, double& x3, double& y3, double& x4, double& y4, double& x5, double& y5, 00173 MsqError& err ); 00174 00175 void basic_tet_test( Vector3D& p3, Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7, Vector3D& p8, 00176 Vector3D& p9, MsqError& err ); 00177 00178 void basic_quad_test( Vector3D& p2, Vector3D& p3, Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7, 00179 MsqError& err ); 00180 00181 void basic_hex_test( Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7, Vector3D& p8, Vector3D& p9, 00182 Vector3D& p10, Vector3D& p11, Vector3D& p12, Vector3D& p13, Vector3D& p14, Vector3D& p15, 00183 Vector3D& p16, Vector3D& p17, Vector3D& p18, Vector3D& p19, MsqError& err ); 00184 }; 00185 00186 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( HigherOrderTest, "HigherOrderTest" ); 00187 CPPUNIT_TEST_SUITE_NAMED_REGISTRATION( HigherOrderTest, "Regression" ); 00188 00189 static inline double dist( double x1, double y1, double x2, double y2 ) 00190 { 00191 return sqrt( ( x1 - x2 ) * ( x1 - x2 ) + ( y1 - y2 ) * ( y1 - y2 ) ); 00192 } 00193 00194 // ideal triangle is defined such that J has unit determinant 00195 const double IDEAL_TRI_SIDE = sqrt( 2.0 / sqrt( 3.0 ) ); 00196 const double IDEAL_TRI_HEIGHT = 1 / IDEAL_TRI_SIDE; 00197 00198 int tri_mid_edge_nodes_edge_center( double x2, double y2, double x3, double y3, double x4, double y4, double x5, 00199 double y5, double epsilon ) 00200 { 00201 double x0 = 0.0, y0 = 0.0; 00202 double x1 = IDEAL_TRI_SIDE, y1 = 0.0; 00203 double x01 = 0.5 * ( x0 + x1 ); 00204 double x12 = 0.5 * ( x1 + x2 ); 00205 double x20 = 0.5 * ( x2 + x0 ); 00206 double y01 = 0.5 * ( y0 + y1 ); 00207 double y12 = 0.5 * ( y1 + y2 ); 00208 double y20 = 0.5 * ( y2 + y0 ); 00209 00210 int result = 0; 00211 if( dist( x3, y3, x01, y01 ) > epsilon ) result |= 1; 00212 if( dist( x4, y4, x12, y12 ) > epsilon ) result |= 2; 00213 if( dist( x5, y5, x20, y20 ) > epsilon ) result |= 4; 00214 00215 return result; 00216 } 00217 00218 void HigherOrderTest::test_tri_basic_ideal() 00219 { 00220 MsqPrintError err( cerr ); 00221 const double eps = 1e-4; 00222 const double x2eq = 0.5 * IDEAL_TRI_SIDE; 00223 const double y2eq = IDEAL_TRI_HEIGHT; 00224 00225 // try starting with the optimal result 00226 double x2 = x2eq; 00227 double y2 = y2eq; 00228 double x3 = 0.5 * IDEAL_TRI_SIDE; 00229 double y3 = 0.0; 00230 double x4 = 0.75 * IDEAL_TRI_SIDE; 00231 double y4 = 0.5 * IDEAL_TRI_HEIGHT; 00232 double x5 = 0.25 * IDEAL_TRI_SIDE; 00233 double y5 = 0.5 * IDEAL_TRI_HEIGHT; 00234 // crit.write_mesh_steps( "test_tri_basic_ideal", TerminationCriterion::VTK ); 00235 basic_tri_test( x2, y2, x3, y3, x4, y4, x5, y5, err ); 00236 // crit.write_mesh_steps( "", TerminationCriterion::NOTYPE ); 00237 ASSERT_NO_ERROR( err ); 00238 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00239 CPPUNIT_ASSERT_DOUBLES_EQUAL( x2eq, x2, eps ); 00240 CPPUNIT_ASSERT_DOUBLES_EQUAL( y2eq, y2, eps ); 00241 int midok = tri_mid_edge_nodes_edge_center( x2, y2, x3, y3, x4, y4, x5, y5, eps ); 00242 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00243 } 00244 00245 void HigherOrderTest::test_tri_basic_mid_spin() 00246 { 00247 MsqPrintError err( cerr ); 00248 const double eps = 1e-4; 00249 const double x2eq = 0.5 * IDEAL_TRI_SIDE; 00250 const double y2eq = IDEAL_TRI_HEIGHT; 00251 00252 // try moving the mid-edge nodes along the edge away from the edge center 00253 double x0 = 0.0, y0 = 0.0; 00254 double x1 = IDEAL_TRI_SIDE, y1 = 0.0; 00255 double x2 = x2eq; 00256 double y2 = y2eq; 00257 const double f = 0.7; 00258 double x3 = f * x0 + ( 1 - f ) * x1; 00259 double y3 = f * y0 + ( 1 - f ) * y1; 00260 double x4 = f * x1 + ( 1 - f ) * x2; 00261 double y4 = f * y1 + ( 1 - f ) * y2; 00262 double x5 = f * x2 + ( 1 - f ) * x0; 00263 double y5 = f * y2 + ( 1 - f ) * y0; 00264 // crit.write_mesh_steps( "test_tri_basic_mid_spin", TerminationCriterion::VTK ); 00265 basic_tri_test( x2, y2, x3, y3, x4, y4, x5, y5, err ); 00266 // crit.write_mesh_steps( "", TerminationCriterion::NOTYPE ); 00267 ASSERT_NO_ERROR( err ); 00268 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00269 CPPUNIT_ASSERT_DOUBLES_EQUAL( x2eq, x2, eps ); 00270 CPPUNIT_ASSERT_DOUBLES_EQUAL( y2eq, y2, eps ); 00271 int midok = tri_mid_edge_nodes_edge_center( x2, y2, x3, y3, x4, y4, x5, y5, eps ); 00272 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00273 } 00274 00275 void HigherOrderTest::test_tri_basic_mid_convex() 00276 { 00277 MsqPrintError err( cerr ); 00278 const double eps = 1e-4; 00279 const double x2eq = 0.5 * IDEAL_TRI_SIDE; 00280 const double y2eq = IDEAL_TRI_HEIGHT; 00281 00282 // try equilateral corners with all egdes convex 00283 double x2 = x2eq; 00284 double y2 = y2eq; 00285 double x3 = 0.5 * IDEAL_TRI_SIDE; 00286 double y3 = -0.3; 00287 double x4 = IDEAL_TRI_SIDE; 00288 double y4 = 0.5; 00289 double x5 = 0.0; 00290 double y5 = 0.5; 00291 basic_tri_test( x2, y2, x3, y3, x4, y4, x5, y5, err ); 00292 ASSERT_NO_ERROR( err ); 00293 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00294 CPPUNIT_ASSERT_DOUBLES_EQUAL( x2eq, x2, eps ); 00295 // CPPUNIT_ASSERT_DOUBLES_EQUAL( y2eq, y2, eps ); 00296 int midok = tri_mid_edge_nodes_edge_center( x2, y2, x3, y3, x4, y4, x5, y5, eps ); 00297 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00298 } 00299 00300 void HigherOrderTest::test_tri_basic_peak_up() 00301 { 00302 MsqPrintError err( cerr ); 00303 const double eps = 1e-4; 00304 const double x2eq = 0.5 * IDEAL_TRI_SIDE; 00305 const double y2eq = IDEAL_TRI_HEIGHT; 00306 00307 // try moving the top vertex up, also move mid-edge nodes proportionally 00308 double x2 = x2eq; 00309 double y2 = 2.0 * y2eq; 00310 double x3 = 0.5 * IDEAL_TRI_SIDE; 00311 double y3 = 0.0; 00312 double x4 = 1.5 * x2; 00313 double y4 = 0.5 * y2; 00314 double x5 = 0.5 * x2; 00315 double y5 = 0.5 * y2; 00316 basic_tri_test( x2, y2, x3, y3, x4, y4, x5, y5, err ); 00317 ASSERT_NO_ERROR( err ); 00318 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00319 CPPUNIT_ASSERT_DOUBLES_EQUAL( x2eq, x2, eps ); 00320 CPPUNIT_ASSERT_DOUBLES_EQUAL( y2eq, y2, eps ); 00321 int midok = tri_mid_edge_nodes_edge_center( x2, y2, x3, y3, x4, y4, x5, y5, eps ); 00322 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00323 } 00324 00325 void HigherOrderTest::test_tri_basic_peak_down() 00326 { 00327 MsqPrintError err( cerr ); 00328 const double eps = 1e-4; 00329 const double x2eq = 0.5 * IDEAL_TRI_SIDE; 00330 const double y2eq = IDEAL_TRI_HEIGHT; 00331 00332 // try moving the top vertex down, also move mid-edge nodes proportionally 00333 double x2 = x2eq; 00334 double y2 = 0.5 * y2eq; 00335 double x3 = 0.5 * IDEAL_TRI_SIDE; 00336 double y3 = 0.0; 00337 double x4 = 1.5 * x2; 00338 double y4 = 0.5 * y2; 00339 double x5 = 0.5 * x2; 00340 double y5 = 0.5 * y2; 00341 basic_tri_test( x2, y2, x3, y3, x4, y4, x5, y5, err ); 00342 ASSERT_NO_ERROR( err ); 00343 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00344 CPPUNIT_ASSERT_DOUBLES_EQUAL( x2eq, x2, eps ); 00345 CPPUNIT_ASSERT_DOUBLES_EQUAL( y2eq, y2, eps ); 00346 int midok = tri_mid_edge_nodes_edge_center( x2, y2, x3, y3, x4, y4, x5, y5, eps ); 00347 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00348 } 00349 00350 void HigherOrderTest::test_tri_basic_peak_over() 00351 { 00352 MsqPrintError err( cerr ); 00353 const double eps = 1e-4; 00354 const double x2eq = 0.5 * IDEAL_TRI_SIDE; 00355 const double y2eq = IDEAL_TRI_HEIGHT; 00356 00357 // try moving the top vertex to the right, also move mid-edge nodes proportionally 00358 double x2 = x2eq + 0.5; 00359 double y2 = y2eq; 00360 double x3 = 0.5 * IDEAL_TRI_SIDE; 00361 double y3 = 0.0; 00362 double x4 = 1.5 * x2; 00363 double y4 = 0.5 * y2; 00364 double x5 = 0.5 * x2; 00365 double y5 = 0.5 * y2; 00366 basic_tri_test( x2, y2, x3, y3, x4, y4, x5, y5, err ); 00367 ASSERT_NO_ERROR( err ); 00368 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00369 CPPUNIT_ASSERT_DOUBLES_EQUAL( x2eq, x2, eps ); 00370 CPPUNIT_ASSERT_DOUBLES_EQUAL( y2eq, y2, eps ); 00371 int midok = tri_mid_edge_nodes_edge_center( x2, y2, x3, y3, x4, y4, x5, y5, eps ); 00372 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00373 } 00374 00375 int quad_all_in_xy_plane( const Vector3D& p2, const Vector3D& p3, const Vector3D& p4, const Vector3D& p5, 00376 const Vector3D& p6, const Vector3D& p7, double epsilon ) 00377 { 00378 Vector3D list[] = { p2, p3, p4, p5, p6, p7 }; 00379 const int n = sizeof( list ) / sizeof( list[0] ); 00380 int result = 0; 00381 for( int i = 0; i < n; ++i ) 00382 if( fabs( list[i].z() ) > epsilon ) result |= ( 1 << i ); 00383 return result; 00384 } 00385 00386 int quad_mid_edge_nodes_edge_center( const Vector3D& p2, const Vector3D& p3, const Vector3D& p4, const Vector3D& p5, 00387 const Vector3D& p6, const Vector3D& p7, double epsilon ) 00388 { 00389 const Vector3D p0( 0.0, 0.0, 0.0 ); 00390 const Vector3D p1( QEL, 0.0, 0.0 ); 00391 const Vector3D e0 = 0.5 * ( p0 + p1 ); 00392 const Vector3D e1 = 0.5 * ( p1 + p2 ); 00393 const Vector3D e2 = 0.5 * ( p2 + p3 ); 00394 const Vector3D e3 = 0.5 * ( p3 + p0 ); 00395 00396 int result = 0; 00397 if( ( p4 - e0 ).length() > epsilon ) result |= 1; 00398 if( ( p5 - e1 ).length() > epsilon ) result |= 2; 00399 if( ( p6 - e2 ).length() > epsilon ) result |= 4; 00400 if( ( p7 - e3 ).length() > epsilon ) result |= 8; 00401 00402 return result; 00403 } 00404 00405 static void get_ideal_quad( Vector3D& p2, Vector3D& p3, Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7 ) 00406 { 00407 p2.set( QEL, QEL, 0 ); 00408 p3.set( 0.0, QEL, 0 ); 00409 p4.set( QEL / 2, 0.0, 0 ); 00410 p5.set( QEL, QEL / 2, 0 ); 00411 p6.set( QEL / 2, QEL, 0 ); 00412 p7.set( 0.0, QEL / 2, 0 ); 00413 } 00414 00415 void HigherOrderTest::test_quad_basic_ideal() 00416 { 00417 MsqPrintError err( cerr ); 00418 const double eps = 5e-2; 00419 Vector3D p2, p3, p4, p5, p6, p7; 00420 const Vector3D p0( 0.0, 0.0, 0.0 ); 00421 const Vector3D p1( QEL, 0.0, 0.0 ); 00422 00423 // try starting with the optimal result 00424 get_ideal_quad( p2, p3, p4, p5, p6, p7 ); 00425 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00426 ASSERT_NO_ERROR( err ); 00427 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00428 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00429 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00430 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00431 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00432 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00433 } 00434 00435 void HigherOrderTest::test_quad_basic_mid_spin() 00436 { 00437 MsqPrintError err( cerr ); 00438 const double eps = 5e-2; 00439 Vector3D p2, p3, p4, p5, p6, p7; 00440 const Vector3D p0( 0.0, 0.0, 0.0 ); 00441 const Vector3D p1( QEL, 0.0, 0.0 ); 00442 00443 // try moving the mid-edge nodes along the edge away from the edge center 00444 p2.set( QEL, QEL, 0 ); 00445 p3.set( 0, QEL, 0 ); 00446 double f = 0.4; 00447 p4 = f * p0 + ( 1 - f ) * p1; 00448 p5 = f * p1 + ( 1 - f ) * p2; 00449 p6 = f * p2 + ( 1 - f ) * p3; 00450 p7 = f * p3 + ( 1 - f ) * p0; 00451 // crit.write_mesh_steps( "quad_basic_mid_spin", TerminationCriterion::VTK ); 00452 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00453 // crit.write_mesh_steps( "", TerminationCriterion::NOTYPE ); 00454 ASSERT_NO_ERROR( err ); 00455 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00456 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00457 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00458 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00459 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00460 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00461 } 00462 00463 void HigherOrderTest::test_quad_basic_mid_convex() 00464 { 00465 MsqPrintError err( cerr ); 00466 const double eps = 5e-2; 00467 Vector3D p2, p3, p4, p5, p6, p7; 00468 const Vector3D p0( 0.0, 0.0, 0.0 ); 00469 const Vector3D p1( QEL, 0.0, 0.0 ); 00470 00471 // try square corners with all egdes convex 00472 get_ideal_quad( p2, p3, p4, p5, p6, p7 ); 00473 p4 += Vector3D( 0.0, -0.2 * QEL, 0 ); 00474 p5 += Vector3D( 0.2 * QEL, 0.0, 0 ); 00475 p6 += Vector3D( 0.0, 0.2 * QEL, 0 ); 00476 p7 += Vector3D( -0.2 * QEL, 0.0, 0 ); 00477 // crit.write_mesh_steps( "quad_basic_mid_convex", TerminationCriterion::VTK ); 00478 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00479 // crit.write_mesh_steps( "", TerminationCriterion::NOTYPE ); 00480 ASSERT_NO_ERROR( err ); 00481 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00482 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00483 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00484 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00485 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00486 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00487 } 00488 00489 void HigherOrderTest::test_quad_basic_left_down() 00490 { 00491 MsqPrintError err( cerr ); 00492 const double eps = 5e-2; 00493 Vector3D p2, p3, p4, p5, p6, p7; 00494 const Vector3D p0( 0.0, 0.0, 0.0 ); 00495 const Vector3D p1( QEL, 0.0, 0.0 ); 00496 00497 // try moving the top left vertex down, also move mid-edge nodes proportionally 00498 get_ideal_quad( p2, p3, p4, p5, p6, p7 ); 00499 p3 -= Vector3D( 0.0, 0.0, QEL / 2 ); 00500 p6 = 0.5 * ( p2 + p3 ); 00501 p7 = 0.5 * ( p0 + p3 ); 00502 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00503 ASSERT_NO_ERROR( err ); 00504 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00505 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00506 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00507 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00508 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00509 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00510 } 00511 00512 void HigherOrderTest::test_quad_basic_top_down() 00513 { 00514 MsqPrintError err( cerr ); 00515 const double eps = 5e-2; 00516 Vector3D p2, p3, p4, p5, p6, p7; 00517 const Vector3D p0( 0.0, 0.0, 0.0 ); 00518 const Vector3D p1( QEL, 0.0, 0.0 ); 00519 00520 // try moving the top two vertices down, also move mid-edge nodes proportionally 00521 get_ideal_quad( p2, p3, p4, p5, p6, p7 ); 00522 p2 -= Vector3D( 0.0, QEL / 2, 0.0 ); 00523 p3 -= Vector3D( 0.0, QEL / 2, 0.0 ); 00524 p5 = 0.5 * ( p1 + p2 ); 00525 p6 = 0.5 * ( p2 + p3 ); 00526 p7 = 0.5 * ( p0 + p3 ); 00527 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00528 ASSERT_NO_ERROR( err ); 00529 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00530 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00531 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00532 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00533 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00534 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00535 } 00536 00537 void HigherOrderTest::test_quad_basic_right_up() 00538 { 00539 MsqPrintError err( cerr ); 00540 const double eps = 5e-2; 00541 Vector3D p2, p3, p4, p5, p6, p7; 00542 const Vector3D p0( 0.0, 0.0, 0.0 ); 00543 const Vector3D p1( QEL, 0.0, 0.0 ); 00544 00545 // try moving the top right vertex up, also move mid-edge nodes proportionally 00546 get_ideal_quad( p2, p3, p4, p5, p6, p7 ); 00547 p2 += Vector3D( 0.0, QEL * 2, 0.0 ); 00548 p5 = 0.5 * ( p1 + p2 ); 00549 p6 = 0.5 * ( p2 + p3 ); 00550 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00551 ASSERT_NO_ERROR( err ); 00552 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00553 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00554 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00555 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00556 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00557 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00558 } 00559 00560 void HigherOrderTest::test_quad_basic_left_over() 00561 { 00562 MsqPrintError err( cerr ); 00563 const double eps = 5e-2; 00564 Vector3D p2, p3, p4, p5, p6, p7; 00565 const Vector3D p0( 0.0, 0.0, 0.0 ); 00566 const Vector3D p1( QEL, 0.0, 0.0 ); 00567 00568 // try moving the top left vertex to the right, also move mid-edge nodes proportionally 00569 get_ideal_quad( p2, p3, p4, p5, p6, p7 ); 00570 p3 -= Vector3D( QEL / 2, 0.0, 0.0 ); 00571 p6 = 0.5 * ( p2 + p3 ); 00572 p7 = 0.5 * ( p0 + p3 ); 00573 basic_quad_test( p2, p3, p4, p5, p6, p7, err ); 00574 ASSERT_NO_ERROR( err ); 00575 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00576 CPPUNIT_ASSERT_EQUAL( 0, quad_all_in_xy_plane( p2, p3, p4, p5, p6, p7, eps ) ); 00577 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( QEL, QEL, 0 ), p2, eps ); 00578 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, QEL, 0 ), p3, eps ); 00579 int midok = quad_mid_edge_nodes_edge_center( p2, p3, p4, p5, p6, p7, eps ); 00580 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00581 } 00582 00583 // ideal tet is defined to have unit lambda (cube root of det == 1) 00584 // const double IDEAL_TET_SIDE = sqrt(2.0) * pow(3.0, 1.0/3.0); // this was for unit area 00585 const double IDEAL_TET_SIDE = 2.0 / pow( 33, 1.0 / 6.0 ); 00586 const double IDEAL_TET_BASE = sqrt( 3.0 ) * 0.5 * IDEAL_TET_SIDE; 00587 const double IDEAL_TET_HEIGHT = sqrt( 2.0 / 3.0 ) * IDEAL_TET_SIDE; 00588 00589 int tet_mid_edge_nodes_edge_center( const Vector3D& p3, const Vector3D& p4, const Vector3D& p5, const Vector3D& p6, 00590 const Vector3D& p7, const Vector3D& p8, const Vector3D& p9, double epsilon ) 00591 { 00592 const Vector3D p0( 0.0, 0.0, 0.0 ); 00593 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00594 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00595 const Vector3D e0 = 0.5 * ( p0 + p1 ); 00596 const Vector3D e1 = 0.5 * ( p1 + p2 ); 00597 const Vector3D e2 = 0.5 * ( p2 + p0 ); 00598 const Vector3D e3 = 0.5 * ( p0 + p3 ); 00599 const Vector3D e4 = 0.5 * ( p1 + p3 ); 00600 const Vector3D e5 = 0.5 * ( p2 + p3 ); 00601 00602 int result = 0; 00603 if( ( p4 - e0 ).length() > epsilon ) result |= 1; 00604 if( ( p5 - e1 ).length() > epsilon ) result |= 2; 00605 if( ( p6 - e2 ).length() > epsilon ) result |= 4; 00606 if( ( p7 - e3 ).length() > epsilon ) result |= 8; 00607 if( ( p8 - e4 ).length() > epsilon ) result |= 16; 00608 if( ( p9 - e5 ).length() > epsilon ) result |= 32; 00609 00610 return result; 00611 } 00612 00613 inline static void get_ideal_tet( Vector3D& p3, Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7, Vector3D& p8, 00614 Vector3D& p9 ) 00615 { 00616 const Vector3D p0( 0.0, 0.0, 0.0 ); 00617 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00618 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00619 p3.set( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00620 p4.set( 0.5 * ( p0 + p1 ) ); 00621 p5.set( 0.5 * ( p1 + p2 ) ); 00622 p6.set( 0.5 * ( p2 + p0 ) ); 00623 p7.set( 0.5 * ( p0 + p3 ) ); 00624 p8.set( 0.5 * ( p1 + p3 ) ); 00625 p9.set( 0.5 * ( p2 + p3 ) ); 00626 } 00627 00628 void HigherOrderTest::test_tet_basic_ideal() 00629 { 00630 MsqPrintError err( cerr ); 00631 const double eps = 5e-2; 00632 const Vector3D p3eq( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00633 Vector3D p3, p4, p5, p6, p7, p8, p9; 00634 const Vector3D p0( 0.0, 0.0, 0.0 ); 00635 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00636 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00637 00638 // try starting with the optimal result 00639 get_ideal_tet( p3, p4, p5, p6, p7, p8, p9 ); 00640 basic_tet_test( p3, p4, p5, p6, p7, p8, p9, err ); 00641 ASSERT_NO_ERROR( err ); 00642 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00643 CPPUNIT_ASSERT_VECTORS_EQUAL( p3eq, p3, eps ); 00644 int midok = tet_mid_edge_nodes_edge_center( p3, p4, p5, p6, p7, p8, p9, eps ); 00645 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00646 } 00647 00648 void HigherOrderTest::test_tet_basic_mid_spin() 00649 { 00650 MsqPrintError err( cerr ); 00651 const double eps = 5e-2; 00652 const Vector3D p3eq( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00653 Vector3D p3, p4, p5, p6, p7, p8, p9; 00654 const Vector3D p0( 0.0, 0.0, 0.0 ); 00655 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00656 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00657 00658 // try moving the mid-edge nodes along the edge away from the edge center 00659 p3 = p3eq; 00660 p4 = 0.6 * p0 + 0.4 * p1; 00661 p5 = 0.6 * p1 + 0.4 * p2; 00662 p6 = 0.6 * p2 + 0.4 * p0; 00663 p7 = 0.6 * p0 + 0.4 * p3; 00664 p8 = 0.4 * p1 + 0.6 * p3; 00665 p9 = 0.3 * p2 + 0.7 * p3; 00666 basic_tet_test( p3, p4, p5, p6, p7, p8, p9, err ); 00667 ASSERT_NO_ERROR( err ); 00668 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00669 CPPUNIT_ASSERT_VECTORS_EQUAL( p3eq, p3, eps ); 00670 int midok = tet_mid_edge_nodes_edge_center( p3, p4, p5, p6, p7, p8, p9, eps ); 00671 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00672 } 00673 00674 void HigherOrderTest::test_tet_basic_mid_convex() 00675 { 00676 MsqPrintError err( cerr ); 00677 const double eps = 5e-2; 00678 const Vector3D p3eq( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00679 Vector3D p3, p4, p5, p6, p7, p8, p9; 00680 const Vector3D p0( 0.0, 0.0, 0.0 ); 00681 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00682 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00683 00684 // try equilateral corners with all egdes convex 00685 get_ideal_tet( p3, p4, p5, p6, p7, p8, p9 ); 00686 p4 += Vector3D( 0.0, -0.2, -0.2 ); 00687 p5 += Vector3D( 0.2, 0.2, -0.2 ); 00688 p6 += Vector3D( -0.2, 0.2, -0.2 ); 00689 p7 += Vector3D( -0.2, -0.2, 0.2 ); 00690 p8 += Vector3D( 0.2, -0.2, 0.2 ); 00691 p9 += Vector3D( 0.0, 0.2, 0.2 ); 00692 basic_tet_test( p3, p4, p5, p6, p7, p8, p9, err ); 00693 ASSERT_NO_ERROR( err ); 00694 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00695 CPPUNIT_ASSERT_VECTORS_EQUAL( p3eq, p3, eps ); 00696 int midok = tet_mid_edge_nodes_edge_center( p3, p4, p5, p6, p7, p8, p9, eps ); 00697 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00698 } 00699 00700 void HigherOrderTest::test_tet_basic_apex_down() 00701 { 00702 MsqPrintError err( cerr ); 00703 const double eps = 5e-2; 00704 const Vector3D p3eq( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00705 Vector3D p3, p4, p5, p6, p7, p8, p9; 00706 const Vector3D p0( 0.0, 0.0, 0.0 ); 00707 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00708 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00709 00710 // try moving the top vertex down, also move mid-edge nodes proportionally 00711 get_ideal_tet( p3, p4, p5, p6, p7, p8, p9 ); 00712 p3 -= Vector3D( 0.0, 0.0, 0.5 ); 00713 p7 -= Vector3D( 0.0, 0.0, 0.25 ); 00714 p8 -= Vector3D( 0.0, 0.0, 0.25 ); 00715 p9 -= Vector3D( 0.0, 0.0, 0.25 ); 00716 basic_tet_test( p3, p4, p5, p6, p7, p8, p9, err ); 00717 ASSERT_NO_ERROR( err ); 00718 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00719 CPPUNIT_ASSERT_VECTORS_EQUAL( p3eq, p3, eps ); 00720 int midok = tet_mid_edge_nodes_edge_center( p3, p4, p5, p6, p7, p8, p9, eps ); 00721 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00722 } 00723 00724 void HigherOrderTest::test_tet_basic_apex_up() 00725 { 00726 MsqPrintError err( cerr ); 00727 const double eps = 5e-2; 00728 const Vector3D p3eq( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00729 Vector3D p3, p4, p5, p6, p7, p8, p9; 00730 const Vector3D p0( 0.0, 0.0, 0.0 ); 00731 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00732 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00733 00734 // try moving the top vertex up, also move mid-edge nodes proportionally 00735 get_ideal_tet( p3, p4, p5, p6, p7, p8, p9 ); 00736 p3 += Vector3D( 0.0, 0.0, 3.0 ); 00737 p7 += Vector3D( 0.0, 0.0, 1.5 ); 00738 p8 += Vector3D( 0.0, 0.0, 1.5 ); 00739 p9 += Vector3D( 0.0, 0.0, 1.5 ); 00740 basic_tet_test( p3, p4, p5, p6, p7, p8, p9, err ); 00741 ASSERT_NO_ERROR( err ); 00742 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00743 CPPUNIT_ASSERT_VECTORS_EQUAL( p3eq, p3, eps ); 00744 int midok = tet_mid_edge_nodes_edge_center( p3, p4, p5, p6, p7, p8, p9, eps ); 00745 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00746 } 00747 00748 void HigherOrderTest::test_tet_basic_apex_over() 00749 { 00750 MsqPrintError err( cerr ); 00751 const double eps = 5e-2; 00752 const Vector3D p3eq( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE / 3.0, IDEAL_TET_HEIGHT ); 00753 Vector3D p3, p4, p5, p6, p7, p8, p9; 00754 const Vector3D p0( 0.0, 0.0, 0.0 ); 00755 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00756 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00757 00758 // try moving the top vertex to the right, also move mid-edge nodes proportionally 00759 get_ideal_tet( p3, p4, p5, p6, p7, p8, p9 ); 00760 p3 -= Vector3D( 0.3, 0.0, 0.0 ); 00761 p7 = 0.5 * ( p0 + p3 ); 00762 p8 = 0.5 * ( p1 + p3 ); 00763 p9 = 0.5 * ( p2 + p3 ); 00764 basic_tet_test( p3, p4, p5, p6, p7, p8, p9, err ); 00765 ASSERT_NO_ERROR( err ); 00766 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00767 CPPUNIT_ASSERT_VECTORS_EQUAL( p3eq, p3, eps ); 00768 int midok = tet_mid_edge_nodes_edge_center( p3, p4, p5, p6, p7, p8, p9, eps ); 00769 CPPUNIT_ASSERT_EQUAL( 0, midok ); 00770 } 00771 00772 void HigherOrderTest::test_hex_basic() 00773 { 00774 MsqError err; 00775 const double P = 0.25; // fraction to perturb higher-order nodes: 0.5->optimal 00776 const double Q = 1 - P; 00777 const double EPSILON = 1e-4; 00778 const unsigned long num_vtx = 27; 00779 00780 // Define a 27-node hex with corner nodes fixed and higher-order 00781 // nodes perturbed a bit. Smooth to see if HO nodes are moved 00782 // to center of edge/face/volume. 00783 double coords[3 * num_vtx] = { // bottom corners 00784 0, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 0, 00785 // top corners 00786 0, 0, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 00787 // bottom mid-edge 00788 P, 0, 0, 1, Q, 0, Q, 1, 0, 0, P, 0, 00789 // vertical mid-edge 00790 0, 0, P, 1, 0, P, 1, 1, Q, 0, 1, Q, 00791 // top mid-edge 00792 Q, 0, 1, 1, P, 1, P, 1, 1, 0, Q, 1, 00793 // vertical faces 00794 P, 0, Q, 1, Q, P, Q, 1, P, 0, P, Q, 00795 // bottom and top mid-face 00796 P, Q, 0, Q, P, 1, 00797 // mid-element 00798 P, Q, P 00799 }; 00800 unsigned long conn[num_vtx]; 00801 for( unsigned i = 0; i < num_vtx; i++ ) 00802 conn[i] = i; 00803 int fixed[num_vtx]; 00804 std::fill( fixed, fixed + 8, 1 ); 00805 std::fill( fixed + 8, fixed + num_vtx, 0 ); 00806 const EntityTopology type = HEXAHEDRON; 00807 unsigned long offsets[] = { 0, num_vtx }; 00808 ArrayMesh one_hex( 3, num_vtx, coords, fixed, 1, &type, conn, offsets ); 00809 00810 // smooth 00811 q.run_instructions( &one_hex, err ); 00812 ASSERT_NO_ERROR( err ); 00813 00814 // test that higher-order nodes were moved to the expected locations 00815 // bottom mid-edge 00816 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 0.0, 0 ), Vector3D( coords + 3 * 8 ), EPSILON ); 00817 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1.0, 0.5, 0 ), Vector3D( coords + 3 * 9 ), EPSILON ); 00818 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 1.0, 0 ), Vector3D( coords + 3 * 10 ), EPSILON ); 00819 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.0, 0.5, 0 ), Vector3D( coords + 3 * 11 ), EPSILON ); 00820 // vertical mid-edge 00821 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 0, 0.5 ), Vector3D( coords + 3 * 12 ), EPSILON ); 00822 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1, 0, 0.5 ), Vector3D( coords + 3 * 13 ), EPSILON ); 00823 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1, 1, 0.5 ), Vector3D( coords + 3 * 14 ), EPSILON ); 00824 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0, 1, 0.5 ), Vector3D( coords + 3 * 15 ), EPSILON ); 00825 // top mid-edge 00826 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 0.0, 1 ), Vector3D( coords + 3 * 16 ), EPSILON ); 00827 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1.0, 0.5, 1 ), Vector3D( coords + 3 * 17 ), EPSILON ); 00828 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 1.0, 1 ), Vector3D( coords + 3 * 18 ), EPSILON ); 00829 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.0, 0.5, 1 ), Vector3D( coords + 3 * 19 ), EPSILON ); 00830 // vertical faces 00831 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 0.0, 0.5 ), Vector3D( coords + 3 * 20 ), EPSILON ); 00832 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 1.0, 0.5, 0.5 ), Vector3D( coords + 3 * 21 ), EPSILON ); 00833 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 1.0, 0.5 ), Vector3D( coords + 3 * 22 ), EPSILON ); 00834 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.0, 0.5, 0.5 ), Vector3D( coords + 3 * 23 ), EPSILON ); 00835 // bottom and top mid-face 00836 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 0.5, 0 ), Vector3D( coords + 3 * 24 ), EPSILON ); 00837 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 0.5, 1 ), Vector3D( coords + 3 * 25 ), EPSILON ); 00838 // mid-element 00839 CPPUNIT_ASSERT_VECTORS_EQUAL( Vector3D( 0.5, 0.5, 0.5 ), Vector3D( coords + 3 * 26 ), EPSILON ); 00840 } 00841 00842 /* 00843 Jason, 00844 00845 Here is a very simple test we should run. 00846 00847 A single quadratic triangle. IMR metric with W = equilateral ideal No geometry. 00848 00849 Nodes: 00850 00851 (x,y) Type F/S/F 00852 ----------------------------------------------------------- 00853 (0,0) Corner Fixed 00854 (1,0) Corner Fixed 00855 (x2,y2) Corner Free 00856 (x3,y3) Mid Free 00857 (x4,y4) Mid Free 00858 (x5,y5) Mid Free 00859 -------------------------------------------------------- 00860 00861 Start with sample points at each of the node. 00862 The optimal element should be equilateral. 00863 Try different initial values for x2,y2,x3, etc. 00864 Try only including some of the sample points to see if still get equilateral. 00865 Probably should add this one to the testSuite. 00866 00867 -- Pat 00868 */ 00869 void HigherOrderTest::basic_tri_test( double& x2, double& y2, double& x3, double& y3, double& x4, double& y4, 00870 double& x5, double& y5, MsqError& err ) 00871 { 00872 // Construct Mesh instance 00873 const int DIM = 3; 00874 const int NVTX = 6; 00875 const int NELEM = 1; 00876 double coords[DIM * NVTX] = { 0., 0., 0, IDEAL_TRI_SIDE, 0., 0, x2, y2, 0, x3, y3, 0, x4, y4, 0, x5, y5, 0 }; 00877 const int fixed[NVTX] = { 1, 1, 0, 0, 0, 0 }; 00878 const unsigned long conn[NVTX * NELEM] = { 00879 0, 1, 2, 3, 4, 5, 00880 }; 00881 ArrayMesh mesh( DIM, NVTX, coords, fixed, NELEM, TRIANGLE, conn, false, NVTX ); 00882 PlanarDomain xy( PlanarDomain::XY ); 00883 00884 // Solve 00885 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &xy ); 00886 q.run_instructions( &mesh_and_domain, err );MSQ_ERRRTN( err ); 00887 00888 // Pass back modified coordinates 00889 x2 = coords[2 * DIM + 0]; 00890 y2 = coords[2 * DIM + 1]; 00891 x3 = coords[3 * DIM + 0]; 00892 y3 = coords[3 * DIM + 1]; 00893 x4 = coords[4 * DIM + 0]; 00894 y4 = coords[4 * DIM + 1]; 00895 x5 = coords[5 * DIM + 0]; 00896 y5 = coords[5 * DIM + 1]; 00897 } 00898 00899 void HigherOrderTest::basic_tet_test( Vector3D& p3, Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7, 00900 Vector3D& p8, Vector3D& p9, MsqError& err ) 00901 { 00902 // Construct Mesh instance 00903 const int DIM = 3; 00904 const int NVTX = 10; 00905 const int NELEM = 1; 00906 const Vector3D p0( 0.0, 0.0, 0.0 ); 00907 const Vector3D p1( IDEAL_TET_SIDE, 0.0, 0.0 ); 00908 const Vector3D p2( 0.5 * IDEAL_TET_SIDE, IDEAL_TET_BASE, 0.0 ); 00909 double coords[DIM * NVTX] = { p0.x(), p0.y(), p0.z(), p1.x(), p1.y(), p1.z(), p2.x(), p2.y(), p2.z(), p3.x(), 00910 p3.y(), p3.z(), p4.x(), p4.y(), p4.z(), p5.x(), p5.y(), p5.z(), p6.x(), p6.y(), 00911 p6.z(), p7.x(), p7.y(), p7.z(), p8.x(), p8.y(), p8.z(), p9.x(), p9.y(), p9.z() }; 00912 00913 const int fixed[NVTX] = { 1, 1, 1, 0, 0, 0, 0, 0, 0, 0 }; 00914 const unsigned long conn[NVTX * NELEM] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9 }; 00915 ArrayMesh mesh( DIM, NVTX, coords, fixed, NELEM, TETRAHEDRON, conn, false, NVTX ); 00916 00917 // Solve 00918 q.run_instructions( &mesh, err );MSQ_ERRRTN( err ); 00919 00920 // Pass back modified coordinates 00921 p3.set( coords + 3 * DIM ); 00922 p4.set( coords + 4 * DIM ); 00923 p5.set( coords + 5 * DIM ); 00924 p6.set( coords + 6 * DIM ); 00925 p7.set( coords + 7 * DIM ); 00926 p8.set( coords + 8 * DIM ); 00927 p9.set( coords + 9 * DIM ); 00928 } 00929 00930 void HigherOrderTest::basic_quad_test( Vector3D& p2, Vector3D& p3, Vector3D& p4, Vector3D& p5, Vector3D& p6, 00931 Vector3D& p7, MsqError& err ) 00932 { 00933 // Construct Mesh instance 00934 const int DIM = 3; 00935 const int NVTX = 8; 00936 const int NELEM = 1; 00937 double coords[DIM * NVTX] = { 0.0, 0.0, 0.0, QEL, 0.0, 0.0, p2.x(), p2.y(), 00938 p2.z(), p3.x(), p3.y(), p3.z(), p4.x(), p4.y(), p4.z(), p5.x(), 00939 p5.y(), p5.z(), p6.x(), p6.y(), p6.z(), p7.x(), p7.y(), p7.z() }; 00940 00941 const int fixed[NVTX] = { 1, 1, 0, 0, 0, 0, 0, 0 }; 00942 const unsigned long conn[NVTX * NELEM] = { 0, 1, 2, 3, 4, 5, 6, 7 }; 00943 ArrayMesh mesh( DIM, NVTX, coords, fixed, NELEM, QUADRILATERAL, conn, false, NVTX ); 00944 PlanarDomain xy( PlanarDomain::XY ); 00945 00946 // Solve 00947 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &xy ); 00948 q.run_instructions( &mesh_and_domain, err ); 00949 ASSERT_NO_ERROR( err ); 00950 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00951 00952 // Pass back modified coordinates 00953 p2.set( coords + 2 * DIM ); 00954 p3.set( coords + 3 * DIM ); 00955 p4.set( coords + 4 * DIM ); 00956 p5.set( coords + 5 * DIM ); 00957 p6.set( coords + 6 * DIM ); 00958 p7.set( coords + 7 * DIM ); 00959 } 00960 00961 void HigherOrderTest::basic_hex_test( Vector3D& p4, Vector3D& p5, Vector3D& p6, Vector3D& p7, Vector3D& p8, 00962 Vector3D& p9, Vector3D& p10, Vector3D& p11, Vector3D& p12, Vector3D& p13, 00963 Vector3D& p14, Vector3D& p15, Vector3D& p16, Vector3D& p17, Vector3D& p18, 00964 Vector3D& p19, MsqError& err ) 00965 { 00966 // Construct Mesh instance 00967 const int DIM = 3; 00968 const int NVTX = 20; 00969 const int NELEM = 1; 00970 double coords[DIM * NVTX] = { 0.0, 0.0, 0, 2.0, 0.0, 0, 2.0, 2.0, 0, 00971 0.0, 2.0, 0, p4.x(), p4.y(), p4.z(), p5.x(), p5.y(), p5.z(), 00972 p6.x(), p6.y(), p6.z(), p7.x(), p7.y(), p7.z(), p8.x(), p8.y(), p8.z(), 00973 p9.x(), p9.y(), p9.z(), p10.x(), p10.y(), p10.z(), p11.x(), p11.y(), p11.z(), 00974 p12.x(), p12.y(), p12.z(), p13.x(), p13.y(), p13.z(), p14.x(), p14.y(), p14.z(), 00975 p15.x(), p15.y(), p15.z(), p16.x(), p16.y(), p16.z(), p17.x(), p17.y(), p17.z(), 00976 p18.x(), p18.y(), p18.z(), p19.x(), p19.y(), p19.z() }; 00977 const int fixed[NVTX] = { 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }; 00978 const unsigned long conn[NVTX * NELEM] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19 }; 00979 ArrayMesh mesh( DIM, NVTX, coords, fixed, NELEM, HEXAHEDRON, conn, false, NVTX ); 00980 00981 // Solve 00982 q.run_instructions( &mesh, err );MSQ_ERRRTN( err ); 00983 CPPUNIT_ASSERT( !hit_iteration_limit() ); 00984 00985 // Pass back modified coordinates 00986 p4.set( coords + 4 * DIM ); 00987 p5.set( coords + 5 * DIM ); 00988 p6.set( coords + 6 * DIM ); 00989 p7.set( coords + 7 * DIM ); 00990 p8.set( coords + 8 * DIM ); 00991 p9.set( coords + 9 * DIM ); 00992 p10.set( coords + 10 * DIM ); 00993 p11.set( coords + 11 * DIM ); 00994 p12.set( coords + 12 * DIM ); 00995 p13.set( coords + 13 * DIM ); 00996 p14.set( coords + 14 * DIM ); 00997 p15.set( coords + 15 * DIM ); 00998 p16.set( coords + 16 * DIM ); 00999 p17.set( coords + 17 * DIM ); 01000 p18.set( coords + 18 * DIM ); 01001 p19.set( coords + 19 * DIM ); 01002 } 01003 01004 /* The mesh consists of a single triangle with 3 corner 01005 and 3 mid-side nodes. Two sides of the triangle are 01006 on the boundary, which consists of two curves (maybe 01007 just straight lines), so one side is on each curve. 01008 The corner vertex at the intersection of the two curves is 01009 fixed, the two mid-side nodes on the boundary are 01010 free (with snap-to), one of the other two corner vertices 01011 is fixed, the other free (with snap-to). The remaining 01012 mid-node vertex is free, being in the interior. To close 01013 the domain, I leave up to you... 01014 01015 ----2* 01016 | | 01017 | | 01018 | --5 4 01019 y2 | | 01020 |y5 | 01021 | | | 01022 ----0*------3-----------1----- 01023 |<--x3->| | 01024 |<-------x1-------->| 01025 01026 * Fixed vertices 01027 01028 */ 01029 void HigherOrderTest::test_tri_open_domain( double& x1, double& x3, double& x4, double y2, double y5, double& y4 ) 01030 { 01031 MsqPrintError err( cerr ); 01032 01033 // Validate input 01034 CPPUNIT_ASSERT( x3 < x1 ); 01035 CPPUNIT_ASSERT( y5 < y2 ); 01036 CPPUNIT_ASSERT( x3 > 0.0 ); 01037 CPPUNIT_ASSERT( y5 > 0.0 ); 01038 01039 // Construct Mesh instance 01040 const int DIM = 3; 01041 const int NVTX = 6; 01042 const int NELEM = 1; 01043 double coords[DIM * NVTX] = { 0, 0, 0, x1, 0, 0, 0, y2, 0, x3, 0, 0, x4, y4, 0, 0, y5, 0 }; 01044 const int fixed[NVTX] = { 1, 0, 1, 0, 0, 0 }; 01045 const unsigned long conn[NVTX * NELEM] = { 01046 0, 1, 2, 3, 4, 5, 01047 }; 01048 ArrayMesh mesh( DIM, NVTX, coords, fixed, NELEM, TRIANGLE, conn, false, NVTX ); 01049 01050 // Construct Geometry 01051 LineDomain xaxis( Vector3D( 0, 0, 0 ), Vector3D( 1, 0, 0 ) ); 01052 LineDomain yaxis( Vector3D( 0, 0, 0 ), Vector3D( 0, 1, 0 ) ); 01053 PlanarDomain xyplane( PlanarDomain::XY ); 01054 DomainClassifier::DomainSet sets[] = { &xaxis, &yaxis, &xyplane }; 01055 01056 vector< Mesh::VertexHandle > verts; 01057 vector< Mesh::ElementHandle > elems; 01058 mesh.get_all_vertices( verts, err ); 01059 ASSERT_NO_ERROR( err ); 01060 mesh.get_all_elements( elems, err ); 01061 ASSERT_NO_ERROR( err ); 01062 CPPUNIT_ASSERT_EQUAL( (size_t)1, elems.size() ); 01063 01064 // Associate mesh with geometry 01065 sets[0].vertices.push_back( verts[1] ); 01066 sets[0].vertices.push_back( verts[3] ); 01067 sets[1].vertices.push_back( verts[2] ); 01068 sets[1].vertices.push_back( verts[5] ); 01069 sets[2].elements.push_back( elems[0] ); 01070 sets[2].vertices.push_back( verts[4] ); 01071 DomainClassifier geom; 01072 DomainClassifier::classify_by_handle( geom, &mesh, sets, 3, err ); 01073 ASSERT_NO_ERROR( err ); 01074 01075 // Solve 01076 MeshDomainAssoc mesh_and_domain = MeshDomainAssoc( &mesh, &geom ); 01077 q.run_instructions( &mesh_and_domain, err ); 01078 ASSERT_NO_ERROR( err ); 01079 CPPUNIT_ASSERT( !hit_iteration_limit() ); 01080 01081 // Pass back modified coordinate values 01082 x1 = coords[DIM * 1 + 0]; 01083 x3 = coords[DIM * 3 + 0]; 01084 x4 = coords[DIM * 4 + 0]; 01085 y4 = coords[DIM * 4 + 1]; 01086 } 01087 01088 void HigherOrderTest::test_tri_open_domain() 01089 { 01090 double x1 = 2; 01091 double y2 = 2; 01092 double x3 = 1; 01093 double y5 = 1; 01094 double x4 = 1; 01095 double y4 = 1; 01096 test_tri_open_domain( x1, x3, x4, y2, y5, y4 ); 01097 CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, x1, 1e-3 ); 01098 CPPUNIT_ASSERT_DOUBLES_EQUAL( 2.0, y2, 1e-3 ); 01099 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, x3, 1e-3 ); 01100 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, x4, 1e-3 ); 01101 CPPUNIT_ASSERT_DOUBLES_EQUAL( 1.0, y4, 1e-3 ); 01102 } 01103 01104 /* 01105 Jason, 01106 01107 Here is an interesting test we should run that's slightly more like the SLAC problem. 01108 01109 A single quadratic triangle. IMR metric with W = equilateral ideal No geometry. 01110 01111 Nodes: 01112 01113 (x,y) Type F/S/F 01114 ----------------------------------------------------------- 01115 (0,0) Corner Fixed 01116 (1,0) Corner Fixed 01117 (0.5,10) Corner Free 01118 (0.50,0.49) Mid Fixed 01119 (0.75,0.5) Mid Free 01120 (0.25,0.5) Mid Free 01121 -------------------------------------------------------- 01122 01123 Start with sample points at each of the node. 01124 01125 The Jacobian at the fixed mid-side node is nearly zero 01126 and the goal is to see if this can be improved by moving the 3 free nodes. 01127 01128 The node at (0.5,1.0) should move downward\u2026. Moving it upward makes 01129 the mid-node Jacobian negative, and the IMR barrier should prevent that. 01130 01131 Try including only some of the sample points to see if all are needed. 01132 01133 -- Pat 01134 */ 01135 void HigherOrderTest::test_tri_slac() 01136 { 01137 MsqPrintError err( cerr ); 01138 01139 // Construct Mesh instance 01140 const int DIM = 3; 01141 const int NVTX = 6; 01142 const int NELEM = 1; 01143 double coords[DIM * NVTX] = { 0.00, 0.00, 0, 1.00, 0.00, 0, 0.50, 10.0, 0, 01144 0.50, 0.30, 0, 0.75, 4.50, 0, 0.25, 4.50, 0 }; 01145 const int fixed[NVTX] = { 1, 1, 0, 1, 0, 0 }; 01146 const unsigned long conn[NVTX * NELEM] = { 01147 0, 1, 2, 3, 4, 5, 01148 }; 01149 ArrayMesh mesh( DIM, NVTX, coords, fixed, NELEM, TRIANGLE, conn, false, NVTX ); 01150 PlanarDomain xy( PlanarDomain::XY ); 01151 01152 // Solve 01153 q.run_instructions( &mesh, err ); 01154 ASSERT_NO_ERROR( err ); 01155 CPPUNIT_ASSERT( !hit_iteration_limit() ); 01156 01157 const Vector3D v0( coords + 0 * DIM ); 01158 const Vector3D v1( coords + 1 * DIM ); 01159 const Vector3D v2( coords + 2 * DIM ); 01160 const Vector3D v3( coords + 3 * DIM ); 01161 const Vector3D v4( coords + 4 * DIM ); 01162 const Vector3D v5( coords + 5 * DIM ); 01163 01164 // Expect vertex 2 to have moved downwards signficantly 01165 CPPUNIT_ASSERT( v2.y() < 2.0 ); 01166 CPPUNIT_ASSERT( v2.y() > 0.5 ); 01167 // Expect it to have stayed in the center along the X axis 01168 CPPUNIT_ASSERT_DOUBLES_EQUAL( 0.5, v2.x(), 1e-3 ); 01169 // Expect free mid-edge nodes to be "between" adjacent corners 01170 CPPUNIT_ASSERT( v4.x() > 0.5 ); 01171 CPPUNIT_ASSERT( v4.x() < 1.0 ); 01172 CPPUNIT_ASSERT( v5.x() > 0.0 ); 01173 CPPUNIT_ASSERT( v5.x() < 0.5 ); 01174 // Expect two free mid-edge nodes to be slightly above the 01175 // line between adjacent corners (slightly convex edge) 01176 // because the fixed mid-edge node on the bottom is moved upward 01177 Vector3D m12 = v2 - v1; 01178 Vector3D m02 = v2 - v0; 01179 double t1 = ( m12 % ( v4 - v0 ) ) / ( m12 % m12 ); 01180 double t2 = ( m02 % ( v5 - v0 ) ) / ( m02 % m02 ); 01181 Vector3D l4 = v0 * t1 * m12; // closed point to v4 on line v1,v2 01182 Vector3D l5 = v0 * t2 * m02; // closed point to v5 on line v0,v2 01183 // Check not to far from line 01184 CPPUNIT_ASSERT( ( l4 - v4 ).length() < 0.1 ); 01185 CPPUNIT_ASSERT( ( l5 - v5 ).length() < 0.1 ); 01186 // Check that edges are convex 01187 CPPUNIT_ASSERT( l4.x() < v4.x() ); 01188 CPPUNIT_ASSERT( l5.x() > v5.x() ); 01189 }