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