MOAB: Mesh Oriented datABase  (version 5.3.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, 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 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines